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This work is devoted to the study of a posteriori error estimation and adaptivity 
in parabolic problems with a particular focus on spatial discontinuous Galerkin 
(dG) discretisations. 

We begin by deriving an a posteriori error estimator for a linear non-stationary 
convection-diffusion problem that is discretised with a backward Euler dG method. 
An adaptive algorithm is then proposed to utilise the error estimator. The 
effectiveness of both the error estimator and the proposed algorithm is shown 
through a series of numerical experiments. 

Moving on to nonlinear problems, we investigate the numerical approximation 
of blow-up. To begin this study, we hrst look at the numerical approximation 
of blow-up in nonlinear ODEs through standard time stepping schemes. We 
then derive an a posteriori error estimator for an implicit-explicit (IMEX) dG 
discretisation of a semilinear parabolic PDE with quadratic nonlinearity. An 
adaptive algorithm is proposed that uses the error estimator to approach the 
blow-up time. The adaptive algorithm is then applied in a series of test cases to 
gauge the effectiveness of the error estimator. 

Finally, we consider the adaptive numerical approximation of a nonlinear 
interface problem that is used to model the mass transfer of solutes through 
semi-permiable membranes. An a posteriori error estimator is proposed for the 
IMEX dG discretisation of the model and its effectiveness tested through a series 
of numerical experiments. 
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Chapter 1 
Introduction 


Partial differential equations are key in the modelling of various physical and 
biological phenomena. The solutions to PDEs are usually unavailable through 
analytical means, so numerical methods are employed in order to approximate 
the solution. Furthermore, a number of nonlinear PDF problems exhibit local 
multiscale behaviour such as boundary or interior layers, interfaces or even local 
space-time blow-up. Such local multiscale features require high local resolution 
of the numerical methods employed in their approximation. Hence, the use of 
numerical methods which can automatically detect and resolve such multiscale 
features is of interest. 

Adaptive algorithms that are driven by a posteriori error estimators he at the 
heart of hnite element analysis. For elliptic problems, there are a wide variety 
of different error estimators available [31 Moreover, adaptive algorithms 

for elliptic problems are relatively well understood; at least for simple elliptic 
problems — see [281IH2] for the standard conforming hnite element method and 
[T8l Ell for the interior penalty dG method. 

For linear parabolic problems, there are many error estimators available in the 
literature for popular discretisations (typically a standard time stepping scheme 
paired with a spatial hnite element discretisation). These error estimators are 
usually composed of an initial condition estimator, a space estimator and a time 
estimator. However, generally speaking, it is unclear how to utilise each of these 
individual estimators to drive adaptivity. Furthermore, while mesh change is 
crucial for the efficient numerical approximation of mobile solutions, it is well 
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known that careless mesh rehnement and/or coarsening can lead to destabilisation 
of the hnite element solution [T^ 15^ . While some progress has been made on the 
construction of adaptive algorithms for parabolic problems [30l 00101102103l 011 
EHl EH EOl |96] , none of the algorithms in the literature have been shown to reduce 
the error estimator that they utilise at the correct rate of convergence with respect 
to the average number of degrees of freedom and the total number of time steps. 

This work will investigate adaptive algorithms for spatial discontinuous 
Galerkin discretisations of parabolic problems with a focus on nonlinear problems. 
To achieve this, we shall look at adaptivity in the context of three different 
problems, to be detailed below. 

In Chapter 2, we introduce notation and state some approximation results and 
general theorems that shall be used throughout this work. We also discuss the 
issue of robustness that arises in the a posteriori error estimation of discontinuous 
Galerkin discretisations of stationary convection-diffusion equations. We then 
introduce a robust error bound for the stationary problem, taken from [HTj, that 
shall be used extensively in this work. 

Chapter 3 deals with linear non-st at ionary convection-diffusion equations. In 
particular, we derive an a posteriori error estimator for a backward Euler dG 
discretisation of the problem. The primary challenge in such a derivation is 
the robustness of the error estimator with respect to the diffusion parameter 
£. We address this through the elliptic reconstruction framework of Makridakis 
and Nochetto [HI] paired with a robust error estimator for the stationary problem 
EH together with a robust treatment of the temporal residual. The proposed 
error estimator can be viewed as the analogue of that given in [104] by Verfiirth 
but with a discontinuous Galerkin spatial discretisation instead of a streamline 
upwind Petrov-Galerkin spatial discretisation. As well as CHI. there are other 
error estimators available in the literature for different discretisations of linear 
non-stationary convection-diffusion equations. In particular, we mention [36] 
wherein the authors produce a robust error estimator through a flux reconstruction 
approach and the work of Picasso and Prachittham [89] wherein they develop an 
error estimator for a Crank-Nicolson temporal discretisation; their discretisation 
and error estimator also provide for the use of anisotropic elements. There are 
also a variety of different a posteriori error estimators available for the 






norm, cf. ^51 l6^ 1^ . In addition to the error estimator we also propose an 
adaptive algorithm, based on that given in [SU] , that utilises different parts of the 
error estimator to control space and time adaptivity. The effectiveness of both 
the error estimator and the proposed algorithm is then tested in four numerical 
experiments. It is also worth noting that adaptive algorithms designed specihcally 
for non-stationary convection-diffusion problems are explored in [SUl lUU] . 

In Chapter 4, we investigate the numerical approximation of blow-up in ODEs. 
More specihcally, we derive an a posteriori error estimator for an ODE with 
polynomial nonlinearity that is discretised using standard time stepping schemes. 
The biggest difficulty in the construction of such an error estimator is having 
to deal with a nonlinear error equation — this can be handled through a local 
continuation argument. A continuation argument is a special type of proof by 
contradiction that is often used to prove existence results for nonlinear parabolic 
PDEs; such arguments have been instrumental in the derivation of a posteriori 
bounds for a variety of nonlinear parabolic problems [131 |37l [71]. A posteriori 
error estimators produced by a continuation argument are conditional in the 
sense that they only hold providing that the estimators involved are sufficiently 
small. In order to use the proposed error estimator to approximate the blow-up 
time, we investigate the design of suitable adaptive algorithms. To that end, 
two adaptive algorithms are proposed and then applied in two test cases under 
different time stepping schemes to compare their effectiveness. Although we 
choose to investigate the numerical approximation of blow-up through a posteriori 
error estimation, there are other ways of approaching this problem that have 
been published in the literature. In [66|, the authors prove existence results for 
numerical approximations to a nonlinear ODE with a polynomial growth condition 
provided that the time step lengths are sufficiently small. For the particular case 
of a polynomial nonlinearity, they show that selecting the time step lengths in a 
certain way yields approach to the blow-up time. In |60|, the authors transform an 
ODE with polynomial nonlinearity through an arc length transformation. They 
then use a forward Euler method to approximate the transformed equation and 
they show that their adaptive algorithm, which is based on their transformation 
plus a tolerance controlled ODE integrator, converges towards the blow-up time 
linearly with respect to the total number of time steps. Finally, in [98| the authors 
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approximate a nonlinear ODE using a ^-method along with a temporal rescaling 
of the ODE and they show that their numerical solution has the same asymptotic 
behaviour as the exact solution. 

In Chapter 5, we build on the results of the previous chapter by investigating 
blow-up in semilinear parabolic PDEs. In particular, we study blow-up in nonlinear 
non-stationary convection-diffusion equations that, based on the results of the 
previous chapter, are discretised using an IMEX dG method. An a posteriori 
error estimator is then proposed for this discretisation of the problem. In order to 
produce a viable error estimator for this problem, there are three major obstacles 
that need to be overcome: the lack of symmetry of the problem, the nonlinear 
error equation and the error due to non-conformity. A posteriori error estimation 
for blow-up in symmetric problems has been considered in [75117^ by Kyza and 
Makridakis, however, such results are not easily generalised to non-symmetric 
problems; the lack of symmetry in the problem can, however, be dealt with 
through the Gagliardo-Nirenberg inequality. The nonlinear error equation is dealt 
with in a similar way to in Ghapter 4 — through a continuation argument and 
the error due to non-conformity can be dealt with via localised bounds for the 
non-conforming part of the error |Sl [701 E]. The proposed error estimator is 
utilised to approximate the blow-up time of the problem through an adaptive 
algorithm that is based on those given in Ghapter 3 and Ghapter 4. The adaptive 
algorithm is then applied to some test problems and the results are compared to 
those given in Ghapter 4. It is worth noting that solution prohles close to the 
blow-up time can also be obtained through the rescaling algorithm of Berger and 
Kohn [m lEH] or the MMPDE method [20l [65]. There is also work looking at the 
numerical approximation of blow-up in the nonlinear Schrodinger equation and its 
generalisations 12 EB EOl Ea [M]. Other numerical methods for approximating 
blow-up in a variety of different nonlinear PDEs can be found in [51 |33l [35l [491 [M] . 

In Ghapter 6, we consider an IMEX dG discretisation of a nonlinear interface 
problem that was introduced in [251 EH] ; based on the works [sii[iaEiiDroni[n2, 
to model the mass transfer of solutes through semi-permiable membranes. There 
are a variety of a posteriori error estimators of both residual and recovery type 
in the literature for different discretisations of interface problems. In particular, 
for the conforming hnite element method [22l [23], the discontinuous Galerkin 


10 






method mi and the hnite volume method jEiEa]. Based on these works, and 
the techniques used in previous chapters, we seek to derive a residual-based a 
posteriori error estimator for this discretisation of the model. The effectiveness 
of the error estimator is then tested through a series of numerical experiments 
utilising the adaptive algorithm that was developed in Chapter 3. 

Finally, in Chapter 7 we summarise the results of this work and discuss ways 
in which this work could be extended. 


11 


Chapter 2 
Preliminaries 


2.1 Sobolev spaces 


Let cj C be a bounded Lipschitz domain with boundary du. For 1 < p < +oo, 
we define the U’ norms by 


v\\lp{ui) 



ess sup |u(a;)| 

X^U) 


and the respective spaces by 


for 1 < p < +CX), 
for p = +CX), 


LP{u) := {u I \\u\\lp(uj) < oo} . 

Note that L‘^{uj) is a Hilbert space with an inner product given by 


{u, v)^ := / uv dx. 


When u is the computational domain H (to be defined later) then because both 
the norm and inner product over H are used frequently in this thesis, the 
relevant subscripts are omitted. Given a multi-index a G the weak derivative 
D°‘ of order |q;| is given by 


D" : = 


(9l“l 

dx'^^dx^^' 
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For k eN and 1 < p < C) 0 , the Sobolev space W^'^{uj) is given by 
:= {u e LP{u) I D^u E LP{u), |a| < k] . 

The spaces W^'P{(jj) are eqnipped with the norms 

\\v\\wk,P{uj) ■= I ^ I 

y|a|<fc J 

||n||^yfe,p(^j) := I |D“n| |x,oo(^) 

|a|<fc 

The space hF^’^(ci;) together with the standard inner prodnct is a Hilbert space 
which we shall denote by := Fractional Sobolev spaces 

in particnlar) are also of significant use when trying to make sense of boundary 
values, in the sense of traces, in Sobolev spaces; we refer to [T] for details. 
Whenever boundary values are used in this thesis, they are to be understood 
in the sense of traces. We define the space H}j{uj), which is the prototypical PDF 
solution space, by 

:= {uEH\u) |n|r, = 0}, 

where Tjj is some subset of du with positive one-dimensional Hausdorff measure; 
if Fd = doj, we denote this space by Hq{u)). Finally, we let denote the 

space of all functions u for which is continuous for all multi-indices a with 

|a| < k. 

For T > 0, the spaces L^(0, T; X) (where X is a real Banach space with norm 
II ■ ||x) consist of all measurable functions n : [0,T] —)■ X for which 

for 1 < p < -f-oo, 
for p = -|-oo. 

We also define H\0,T;X) := {u E L^{0,T;X) \ Ut E L2(o,T;X)}. Finally, we 
denote by C(0, T ; X) and C°’^(0, T; X), respectively, the spaces of continuous and 


\v\\lp{0,t-,x) ■= \\v{t)\\P^dt^ < oo 

\v\\lp(0,t-,x) ■= ess sup ||n(t)||x < oo 

0 <t<T 


for 1 < p < -l-cxD, 
for p = -|-cx). 
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Lipschitz continuous functions n : [0,T] —)■ X such that 




max ||n(t)||x < oo, 

0<t<T 


fI |coa(o,T;X) 


:= max 


|'l^||c(0,r;X), 


dv 

m 


< oo. 


L^{0,T-,X) 


2.2 Stationary convection-diffusion equation 

Let the computational domain C be a polygon with boundary this 
assumption will be used throughout the rest of this thesis. We consider the model 
problem of hnding u : —)• M such that 

—eAu + a ■ Vu + bu = f in 

( 2 . 1 ) 

u = 0 on dfl. 


The variable functions are collectively referred to as the data of the problem. 
Problem (2.1) is the prototypical convection-diffusion equation. If the convection 


is constant then the scale of the solution to (2.1) can be characterised through 
the ratio of convection to diffusion as described by the Peclet number 


Pe := 


a||fl| 


When Pe 3> 1, ( 2.1[ ) is advection dominated and can exhibit some or all of the 
following features (see [69l [93] for a detailed analysis): 


• The presence of ordinary layers, spatial areas containing steep gradients of 
the solution u of width 0(e) that usually occur near the outflow boundary 
as boundary layers. 


• The presence of parabolic layers, spatial areas containing moderate gradients 
of the solution u of width O {y/e) that typically occur on the inflow boundary 
as boundary layers or as interior layers. 

These complex spatial features can render the numerical approximation of the 
solution difficult as discussed in the next section. 
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In order to analyse (2.1) and its parabolic counterparts, we must make some 
assumptions on the data. We assume that: 0 < £ < 1, / G a G 

and b G L°°{Q). Furthermore, we require some additional assumptions in order 
to state standard coercivity and continuity results (see, e.g., [971 1110] )- To that 
end, we assume that there are constants (3 >0 and c* > 0 such that 


b - -V ■ SL> (3 a.e. in 11, 116 - V ■ a| \L°°{n) < c^(3- (2.2) 


The weak form of (2.1) reads: hnd u G Hq{Q) such that 


B (u, v) = (/, v) Vv € 


(2.3) 


where 


B(u, v)= {eVu ■ Vv + a ■ Vuv + buv) dx. 

Jn 


2.3 Discontinuous Galerkin method 


(2.4) 


A finite element method is a numerical technique for hnding approximate solutions 
to the weak formulation of PDFs characterised by the use of a subdivision of H 
referred to as the mesh or triangulation. The mesh is a collection of elements 
with K denoting a generic element. A hnite element space is then constructed 
over the mesh and the weak form of the PDF is discretised; different choices of 
hnite element space and different discretisations give rise to different hnite element 
methods. The discretisation parameters are quantities related to convergence of 
the method, specihcally, the diameters of elements in the mesh (and the lengths 
of time steps in parabolic problems). 


When it comes to the hnite element approximation of (2.3), the presence 


of layers introduces a certain amount of difficulty. In particular, the standard 
conforming hnite element method performs poorly if an insufficient number of 
elements are placed in the vicinity of the layers resulting in unphysical oscillations. 
This issue can be solved with layer-adapted meshes such as Shishkin meshes (see 
I7g for an overview of this subject) but these special meshes require a priori 
knowledge of where the layer will occur. These problems led to the development of 
stabilised hnite element methods for convection-dihusion equations [93] the most 
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popular of which are the streamline upwind Petrov-Galerkin (SUPG) method [19] 
and the discontinuous Galerkin (dG) method. In this work, we focus upon a 
dG discretisation of (2.3) taken from [62] which is based upon a classical interior 


penalty discretisation of the diffusive term originally introduced in [71 [9l ETj and 
an upwind discretisation of the transport term first discussed in 


In order to state the dG discretisation of (2.3), we need some additional 
notation. The mesh ( is assumed to be constructed via affine mappings 
Fk : K ^ K with non-singular Jacobian where K is the reference triangle or 
the reference square. The mesh is allowed to contain a uniformly fixed number of 
regular hanging nodes per edge. We define the finite element space 


14 = 14(0 := {u e \v\KoFKeV^{k), iPeC}, (2.5) 


where V^{K) is the space of polynomials of total degree p if iP is the reference 
triangle, or the space of polynomials of degree p in each variable if K is the 
reference square. Let T(0 denote the set of all edges in the mesh C, and 
the set of all interior edges. We also denote the diameter of an element K E ( 
by hx and the length of an edge E E S{() by hx- The outward unit normal to 
the boundary of an element K is denoted by nx- We assume that the mesh ( is 
shape-regular, that is, there exists G > 0 such that for all iL G C we have 


where dx denotes the diameter of the largest ball that can be completely contained 
in K. 

In what follows, it will be useful to associate patches with each element K E (. 
Specifically, we have the (elemental) patch K which is the union of all elements 
that “neighbour” K and the edge patch Ke which is the union of all edges that 
intersect the boundary of K. Formally, these are defined, respectively, by 


K' EC 

kE -.= {[jE, eeS{c) 


dKndK' ^ 0 |, 
aiFGG ^ 0|. 


We also associate an (elemental) patch E with each edge E E £{() which is the 
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union of all elements whose boundary intersects E, viz., 


E:={\Jk, KeC 


EndK^ 


Given an edge E G shared by two elements K and K', a vector held 

V G and a scalar held v G we dehne jumps [■] and averages 

{■} of V and v across E by 


'\k+^\k'): 

[v] := v4 ■ rii^ + v4, ■ 

\k + 

[u] := + vIkiUk'. 


2 

1 


If i? C dfl, we set {v} := v, [v] := v ■ n, {u} := v and [u] := un, with n denoting 
the outward unit normal to the boundary dQ. 

We dehne the inhow and outhow parts of the boundary dfl, respectively, by 

dflin := {x G I a(x) ■ n(x) < 0}, dQout ■= {x G 9G | a(x) ■ n(x) > 0}. 

Similarly, the inhow and outhow parts of an element K are dehned as 

dKin := {x G dK \ a(x) ■ rii^(x) < 0}, dKout '■= ^ dK \ a(x) ■ n^ix) > 0}. 


With all the above notation at hand, the dG approximation to (2.3) reads as 
follows: hnd u/i G 14 such that 


B{uk, lift) + h\{uk, VI,) = (/, lift) Viift e 14, 


( 2 . 6 ) 


where 


B{uh, Vh) := X' / - SiUh) ■ Vvh + (fe - V ■ ai)uhVh dx 

7£ 

Je 


+WftWds, ^ 2 . 7 ) 


E&S{0 


Kec 


'dKout 


Kh{uh,Vh) ^ / {eVuh} ■ [vt] + • [uh] ds. 

EeSiO 
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The penalty parameter, 7 , is set to 7 = 2p‘^ in light of [62] so that the operator 
B + Kh is coercive on I 4 (see below). 

We note that the bilinear form is not well-dehned for arguments in HliVl), 


but the bilinear form B is and is equal to that appearing in (2.3). To analyse the 
dG discretisation, we introduce the quantities 


mm := I (^ll'^“llL 2 (i^) + - 


m 


LHK) 


)+ ||[m]||^ 2 (£;) 


1/2 


E &£{0 


u A : = 


sup 
veH^{n)\{o} 


Lau-Vvdx\ ^ hB..r ...2 

—1 + E viimiiip(Ei 


1/2 


Ee£(0 


These quantities dehne norms on hfo(G) + I 4 . In the literature, ||| • ||| is referred 
to as the energy norm while the quantity | ■ |a is referred to as a dual norm. It is 
easy to see that the bilinear form B is coercive on Hq{Q), viz.. 


B{y,v) > 11 |n| 


(2.8) 


for all V G HliVl), and is continuous in the following sense 


B{u,v) < (|||m||| + |KU)|||n|||, (2.9) 

for all u G Hq{Q) + I 4 and v G Hq{Q). Moreover, the discrete bilinear form is 
coercive for Vh G 14 with respect to the energy norm, viz., 

B{vh,Vh) + Kh{vh,Vh) > ||4ft|ir (2.10) 


The symbols < and > used above and throughout the rest of the thesis are used 
to describe inequalities that are true up to an unspecihed positive constant that 
is independent of the data, the discretisation parameters, the exact solution and 
the dG solution. 
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2.4 Error bounds for the stationary problem 

Let u be the exact solution of a PDE and Uh be some finite element approximation; 
an a posteriori error estimator, rj, is an approximation of the error, e := u — Uh, in 
a certain norm || ■ || such that ||e|| ~ rj. The error estimator must be computable 
and thus is allowed to depend upon the data, the discretisation parameters and the 
hnite element solution Uh but not the unknown solution u. In order to discuss how 
good T] is at approximating ||e||, it is useful to introduce the notion of reliability 
and efficiency of an estimator. An a posteriori estimator, rj, is said to be reliable 
if there is C > 0, independent of the exact solution u, such that 


<Cv, 


( 2 . 11 ) 


while 7] is said to be efficient if there is c > 0, independent of the exact solution 
u, such that 

CT] < \\e\\. (2-12) 


The constants appearing in (2.11 )-(2.12) are often impossible to calculate explicitly 


which leads us to the useful notion of the effectivity index. If u is known for 
particular data, ||e|| may be calculated explicitly for different realisations of Uh- 
Thus, we can compute the effectivity index - the ratio of rj to ||e||: 


effectivity index : = 


V 


The effectivity index naturally leads to the notion of robustness. An error 
estimator, rj, is said to be robust (with respect to || • ||) if the constants in 


(2.11)-(2.12) are always independent of the data, the discretisation parameters 


and the finite element solution Uh. There is also the weaker notion of asymptotic 
robustness: rj is said to be asymptotically robust if it is robust once the discretisation 
parameters are sufficiently small. If rj is asymptotically robust then it successfully 
reproduces the convergence rate of | |e| | with respect to the discretisation parameters. 


Regarding the a posteriori error estimation of (2.3), many estimators exist for 


stabilised hnite element schemes O El El HHl IMl EH El EH El 11081111011113] and 
the primary issue is robustness with respect to the small parameter £. Before the 
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Figure 2.1: Numerical approximation of a boundary layer in the pre-asymptotic 
regime (left) and the asymptotic regime (right). 

layers have been covered by a sufficient number of elements, the pre-asymptotic 
regime, the error in the energy norm is relatively constant. When a sufficient 
number of elements have been placed within the layers, the asymptotic regime, the 
error in the energy norm starts to display the correct convergence rate with respect 
to the discretisation parameters. In order for an error estimator to be robust with 
respect to £ in the energy norm, it must accurately capture the behaviour of 
the energy norm in both the pre-asymptotic and asymptotic regimes. Standard, 
classical estimators in the literature signihcantly overestimate the error in the 
energy norm in the pre-asymptotic regime to allow the layers to be detected and 
rehned but this means that they are not robust with respect to £ in the energy 
norm until the layers have been sufficiently resolved. One way out of this difficulty 
is to add an additional term, the dual norm, to the energy norm to account for 
the pre-asymptotic regime [SI [95] ; robustness is then recovered in all regimes for 
the full norm [971 1110] . The question of whether a robust error estimator exists 
for the energy norm in all regimes is still open, however, recent results in this 
direction seem promising [53] ■ 

An a posteriori estimator for the stationary problem, inspired by m, will 
be utilised in our analysis. More specifically, we have the following result whose 
proof is completely analogous to Theorem 3.2 in [97] and is therefore omitted for 
brevity. 
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Theorem 2.1. For f G let G Hq{Q,) be such that 


B{u% v) = (/, v) Mv G Hl{Q), 


and consider G Vh such that 


B{ul, Vh) + Kh{ul, Vh) = (/, Vh) 'ivh G Vh- 


Then the following a posteriori bound holds for any 0 7 ^ u G Hq^Q): 


B{u^ — ul, v) 






E&SiO 


1/2 


h 






£;e£*"*(C) 


2.5 Finite element approximation results 

Error bounds for the approximation of functions in Hf,{Q) have been constructed 
in the hnite element literature and will be utilised below. 

Theorem 2.2. Given u G Hj^iQ,), there is a finite element quasi-interpolant 
Ixu G Hf,{'[l) n Vh such that 

h~K\\u- IxuWl^^k) < \ \^u\\^2^K) WK G C, 

hE^\\u - IxuWl^e) < l|VM||i2(E) VE G ^(C). 

Proof. See |109] . □ 

A useful tool in a posteriori error estimation involving dG hnite element spaces 
is the approximation of functions in the dG hnite element space, 14, by functions 
in the conforming hnite element space, He{Q) fl I 4 . The next theorem gives 
bounds on such approximation errors. 

Theorem 2.3. Given Uh G I 4 , there exists a decomposition Uh = Uh^c + Uh,d with 
Uh,c £ He{Q) n 14 and Uh,d £ 14 such that the following bounds hold for each 
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element K G (: 


E '%'IINIIi>(E,+ E 

EcKE\dn EcKEnTo 

\\'^h,d\\'L2(^K) < E E hE\\Uh\\‘L2{E)y 

EcKE\dn EcKEnTo 

\\Uh,d\\L°°{K) ^ ll['^/i]llL°°(A'B\an) + \ \'^h\\L^(^Kj^r)To)- 

Proof. The proof is based upon a weighted averaging of the dG solution around 
the nodes; see [zniin] for the first two estimates and [31] for the final estimate. □ 

Given a function u G the projection of u onto 14, denoted by 

IhU G 14 , is the unique solution of 


(n - hu, Vh) = 0 Vuft G 14 . 


(2.13) 


Remark 2.1. For the dG finite element space, IhU can also be constructed by 


considering (2.13) elementwise instead of globally and then piecing together the 
resultant local functions. 


The final theorem in this section gives bounds on the approximation error 
involved in the projection. 

Theorem 2.4. Given a function u G and its projection IhU G 14, the 

following bound holds for any element K E (: 


hK\\u - IhuWi^K) < ||Vn| 42 (^). 

Proof. We use the definition of the projection along with the Gauchy-Schwarz 
inequality to conclude that for any G 14: 

I I'll Ih ^11 L'^i^K) dihU, 'll IfiU^K 

= {U- IhU,u)K 

= {U- IhU,U- Vh)K 

< \\u — IhU\\L'^{K)\\u — Vh\\L'^(K)- 
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Therefore, 


\\u — IhU\\i,2i^K) < \ \u — Vh\\L'^(K)- 

All that remains is to choose Vh to give the bound presented in the theorem. In 


particular, the hnite element interpolant in Theorem 2.2 suffices. 


□ 


2.6 Useful inequalities 

In this section, we introduce general theorems that will be used throughout the 
rest of this thesis. 


Theorem 2.5 (Young’s inequality). Given a, 6 G M, then for any 6 > 0 we 

have 


ah < 


a^5 62 

~Y^Y5' 


Proof. Follows by expanding the inequality — b6 > g_ 


□ 


Theorem 2.6 (Power mean inequality). Given a, 6 > 0, then for any /i > 0 

we have 

(a + 6)'' <max{l,2'^-^}(a^ + 6^). 

Proof. Follows from Jensen’s inequality. □ 


Theorem 2.7 (Multidimensional integration by parts). Given u G 


and V G H{div, fl) := < v G [//^(q)] 


V • V G 


we have 


V ■ Vu dx+ mV ■ v dx = 


MV 


nds. 


I an 


Proof. Follows from the divergence theorem. 


□ 


Theorem 2.8 (Poincare-Friedrichs inequality). Foru G P[f){Q), we have the 
following bound: 

ll«|p <I|Vm||2. 

Proof See |112] . □ 
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Theorem 2.9 (Gagliardo-Nirenberg inequality). For any u G and 

> 0, we have the following bound: 

II'“IIL2+M(r2) ~ ii«iniv«ir 

Proof. See [2] specifically or [HH] for a larger class of inequalities. □ 

Theorem 2.10 (Trace inequality). Givenu G the following bound holds 

for any 6 G (0,1).- 


II“IIl2(9o) ~ '5||V'u||^ + (i + (5 iimII^. 

Proof. See Theorem 1.5.1.10. in [58]. □ 

Theorem 2.11 (Inverse estimate). Given Vh G 14, the following bound holds 
for any K E C- 

Proof See [59] . □ 

The final two theorems in this section relate to a specific group of integral 
inequalities and provide an upper bound on the functions that satisfy such 
inequalities. The first theorem is (classical) Gronwall’s inequality while the second 
theorem is a variant of the first theorem that makes use of a lower order term. 

Theorem 2.12 (Gronwall’s inequality). Let T > 0 and suppose that 
Co, Cl G L^(0,T) and u G 1T^’^(0,T). If for almost every t G (0,T] we have 

u'{t) < Co{t) + ci{t)u{t), 


then 


where G(s, t) 
a.e. then 


m(T) < G(0,T)m( 0) + / G(s,T)co(s)ds, 


exp I / ci(.^)(i.^ . Additionally, if cq and ci are non-negative 


M(r) < G(0 , T) ( m(0 ) + J co(s)(is). 
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Proof. Multiplying the inequality by G ^0,^) yields 


G ^(0, t)M'(t) < G ^(0, t)co(t) + G ^{0,t)ci{t)u{t). 


The product rule and the fundamental theorem of calculus imply that 
^ (G“^(0,f)u(f)) = G~^{0,t)u'{t) — G~^{0,t)ci{t)u{t). 

(JjL 

Thus, 

^{G-\0,t)u{t)) < G-\0,t)co{t). 

The primary result then follows by integrating over [0,T] and noting that 
G(s,T) = G(0,T)G-H0,s). □ 


Theorem 2.13. Let T > 0 and suppose that Cq is a constant, Ci and C 2 are 
non-negative functions and that u is a non-negative function that satisfies 

u^(T) < Cq + f ci{s)u{s) ds-\- f C 2 {s)u‘^{s) ds, 


then 


u{T) < ( |co| + V / ci(s) ds ) exp ( - / C 2 (s) ds ). 


Proof. See Theorem 21 in [35] . 


□ 
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Chapter 3 

A posteriori error estimation and 
adaptivity for non-stationary 
convection-diffusion problems 


3.1 Non-stationary convection-diffusion equation 

For T > 0, we consider the model problem of finding m : x (0, T] —)■ M snch that 


du ^ ^ , r 

— - ei\u + a ■ Vm + bu = j 

ot 

u = 0 

m(-,0) = Uq 


in X (0,T], 

on dfl X (0, T], 
in ff. 


(3.1) 


The model problem (3.1) can display the same spatial featnres as (2.1), namely, 


bonndary and interior layers; the addition of the time domain has the potential 
to make things more complicated, however. Indeed, for snitable data the solution 
may exhibit moving boundary or interior layers. Additionally, the temporal 


behaviour of (3.1) can also be strongly influenced by £ 


The standard weak formulation (see, e.g.. 


of (3.1) reads: hnd u G 


L'^(0,T; n if^(0,T; such that for almost every t G (0,T] we have 

= {f,v) Vn G ilo(fl), (3.2) 
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where B(t]-^-) denotes the bilinear form (2.7) with the data evaluated at the 
time t. We also make some assumptions on the data: Uq G Hq{Q), 0 < e < 1, 
a G [^(0, T; b G ^(0, T; and / G ^(0, T; L\Q)). 

In order to ensure coercivity and continuity of for any t G [0,T], we 


must extend (2.2). To that end, we assume that there are constants f3 > 0 and 


c* > 0 such that 


b--V-SL>/3 a.e. in n X [0,T], ||&- V ■ a||c(o,T;L°°(n)) < c*/?. (3.3) 


3.2 Space-time discretisation 


The semi-discrete discontinuous Galerkin approximation to (3.2) then reads as 


follows. For t = 0, set Uh{0) G 14 to be some projection of uq onto 14- Then, seek 
Uh G 4°’^(0,T; 14) such that for almost every t G (0,T] we have 


duh 

dt 


,Vh] + B{t-,Uh,Vh) + Kh{uh,Vh) = {f,Vh) Vufe G 14- (3.4) 


We shall also consider a full discretisation of problem (3.2) by using a backward 
Euler method to approximate the time derivative. 


To this end, consider a subdivision of [0,T] into time intervals of lengths 


7-1, 




such that = T for some n > 1 then set := 0 and := ■ 


i=i 


i=i 

-k 


Denote an initial triangulation by and then associate a triangulation C,^ to 
each time step fc > 0 which is assumed to have been obtained from by 
locally rehning and coarsening This restriction upon mesh change is made 

to avoid altering the mesh too much between time steps in an attempt to prevent 
degradation of the hnite element solution, cf. [T^lEn]- To each mesh we assign 
the hnite element space := 14(C^) given by (2.5). We also set := /(.,t^). 


:= a(.,t^) and b^ := b[.,t^) for brevity. 

The fully-discrete dG method then reads as follows. Set to be a projection 
of Uo onto V^. For k = 0, ..., n — 1, hnd G such that 


H _(4 


„.A;+1 


+ ') = (3.6) 
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for all G We shall take u\ to be the orthogonal projection of uq 

onto 14 ° although other projections onto can also be used. 


3.3 Error bounds for the non-st at ionary problem 


Here, we will present a posteriori error bounds for both the semi-discrete and 
fully-discrete schemes which can be found in [27]; these results are extended by 
the authors in [72] to convection-diffusion problems with nonlinear reaction term. 


Other a posteriori error estimators for different space-time discretisations of (3.2) 
can be found in [Ml EH] ES] EH] IQOHTTH]. 

To devise our error bounds, we use the elliptic reconstruction approach originally 
introduced by Makridakis and Nochetto for the conforming hnite element method 
[8T] and extended to dG methods in [55l EE] ; this effectively allows decomposition 
of the error into separate parabolic and elliptic parts and can be viewed as 
the a posteriori counterpart to Wheeler’s elliptic projection jlllj from a priori 
analysis. Elliptic reconstruction allows us to bound the elliptic part of the error 
with any error estimator for the stationary problem (2.3) that currently exists in 
the literature; we shall use the bound in Theorem 2.1 taken from [97]. Given that 
we already have a reasonable spatial estimator to use, it is clear that the main 
challenge in the a posteriori estimation of (3.2) is thus related to the parabolic part 
of the error. In particular, the primary challenge is obtaining an error estimator 
that is either robust or at least asymptotically robust with respect to e in the 
type norm 

1/2 


u * := 


\u\ 


L°°(0,r;L2(n)) 


+ 


\u\ 


dt 


This requires careful treatment of the temporal part of the residual and we 
introduce a novel way of dealing with this difficulty. 


3.3.1 An a posteriori bound for the semi-discrete method 

To highlight the main ideas, we begin by deriving an a posteriori bound for the 
semi-discrete method. 
















Definition 3.1. For almost every t G (0,T], we define the elliptic reconstruction 
w G Hq{Q) to be the (unique) solution of the problem 

B{t;w,v) = \/veH({n). 

Remark 3.1. The dG discretisation of the above equation is to find a function 
Wh G T, i4) such that for almost every t G (0,T] we have 


B{t; Wh, Vh) + Kh{wh, Vh) = { f 


duh 

dt 


■,Vh 


yvh G Vh. 


In conjunction with (2.10) and (3.4), this implies that Wh = Uh- Therefore, 


B{t-,w — Uh,v) can be estimated using Theorem 2.1 


Theorem 3.1. The dG solution, Uh, admits a decomposition into a conforming 
part Uh,c £ Hlifil) n Vh and a non-conforming part Uh,d G Vh with Uh = Uh,c + Uh,d 
such that 

(Ill-will+ I«mU)"< E ()(^ + /3'>e)|IWII1(e)+ E 


du 


h,d 


dt 


E 6 £(C) 

^ ^ hE 

E6£(C) 


E 6 £-(C) 


duh 

dt 


L2(E) 


Proof. The first estimate follows from the definition of the norms, Theorem |2.3 
and the Cauchy-Schwarz inequality, viz.. 


(|||Mh,d||| + \Uh,d\A)‘^ < ^ + fi\\Uh,d\\L'2{K) + ^ 


KeC 


E 6 £(C) 


E hr + /9ftdllMllhE)+ E TllWlllw 


hf 


hr 


E&£{0 


< E 

E&£iO 




E&SiO 


The second estimate follows directly from Theorem 2.3 


□ 


The error, e := u — Uh, is decomposed into a parabolic part p and an elliptic 
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part e, viz., 


e = p + e, with p := u — w and e := w — Uh- 


We further dehne Cc ■= u — upc and ec := w — upc to help facilitate the error 
analysis. We are now ready to state our a posteriori error estimator, p, given by 


p : = 


|e(0) I r + ^ df + min I ps2 dt^ , «t f Vs2 dt } + max 


'0 


0<t<T 


1/2 


where ax ■= min with 


'(I. - E f 


Kei 


_ ^ + sAuh - a ■ Vuh - buh 
at 




E&£{0 




hf 


hr 


E6£:(C) 


hi - 


Y1 

E&SiO 


duh 

dt 


LHE) 


hi •“ dE\\[Uh\\\‘i2(^E)- 

EesiC) 


Theorem 3.2. The error the semi-discrete dG method (3.4) satisfies the bound 


lelL < 


p. 


Proof. We know that the exact solution satishes 


(^w) +5(f;M,'t;) = (Z,!;) E H^{n). (3.6) 


Using Dehnition 3.1 we obtain 


+S(f;p,n) = 0 \/v E H^{n). 


(3.7) 
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Setting n = Cc in the above equation gives 


ec j + B(t; Cc, ej = ( ) + B(t; e, ej + B(t; Uh.d, Cc). (3.8) 


Using the elliptic reconstruction property of e along with Theorem 2.1 yields 


B(t;e,e^) < i)s,llkt|||, 


(3.9) 


while using the continuity of the bilinear form B with Theorem 3.1 gives 


B(t;uh,dyec) < (llluft.rflll + k/,,dU)|||ec||| < r/sillkc 


(3.10) 


Combining these results and using the coercivity of the bilinear form B, the 


Cauchy-Schwarz inequality, Young’s inequality and Theorem 3.1 gives 


d 


^(Ikcir^) + IlkcllP 5 4, + >)Salkk|. 


(3.11) 


Let To G [0,T] be such that Ec := ||ec(To)|| = \\ec\\L°°(o,T-,L‘^(n)) then integrating 


(3.11) over [0,To] and using Young’s inequality gives 


Eci\\ec(0)ff+ rflit+l / 173 , * 


(3.12) 


Going back to (3.11) and integrating over [0,T] yields 


le^IlP (it < ||e^(0)|P + / r]l^dt + Ec[ [ rjs^dty 


(3.13) 


Adding (3.12) and (3.13) then using Young’s inequality we obtain 


< 


“ + / A / Os, dt I . 


(3.14) 


Going back to (3.11), integrating over [0,To] and [0,T] then summing the results 
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yields 


kcllj <l|ec(0)|P+ f jfs.dt+ f i(sj|e„||(i«. 


(3.15) 


^0 ^0 

Observing that ||ec|| is contained in |||ec||| up to the coefficient f3, then using the 
Cauchy-Schwarz inequality and Young’s inequality gives 


^C\ I* 


r-T 

+ / Vk dt + I vk dt. 

Jo Jo 


(3.16) 


Going back to (3.15), we can instead use the PoincarGFriedrichs inequality to 


bound I led I by ||Ved| then observe that ||Ved| is also part of |||ed|| up to the 
coefficient e, thus using the Cauchy-Schwarz inequality and Young’s inequality 


gives 


|2 < 


+ [ vkdt + e ^ [ ril^dt. 
Jo Jo 


Putting (3.14), (3.16) and (3.17) together we obtain 


(3.17) 


< 


+ / hli dt + min 
Jo ^ 


VS2 dt 


.at / 


(3.18) 


Obviously from the triangle inequality we have 


P < 




(3.19) 


Thus, all we need to do to complete the proof is to bound and ||ec(0)|p. 

To bound ||'U/i,d||^, we use the dehnition of the norm along with Theorem 2.3 and 
Theorem |3.1[ viz.. 


I kh,o 


pT pT 

Y ' \\\uh,d\\\‘^ dt + max \\uh,d\\‘^ < ‘ ^ - - „2 


0<t<T 


r]Si dt + max r]i. (3.20) 


0<t<T 


Finally, ||ec(0)|p is bounded using the triangle inequality and Theorem 2.3 


VIZ., 


2 < 


e(0)||^ max \\uh,d\\^ < ||e(0)||^ max r/l. (3.21) 


0<t<T 


0<t<T 


These bounds complete the proof. 


□ 
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3.3.2 An a posteriori bound for the fully-discrete method 


We now continue by applying to the fully-discrete setting the general framework 
presented in the previous subsection. 


Definition 3.2. We define G 14^ to he the unique solution of the elliptic 
problem 


B{t^ 


ulwl) +Kh{u 


h-: '^h 






Vt 


Remark 3.2. Fork > 1 we obtain from (3.5) that 
where is the Lf projection operator onto . 


u 


k+l 






Definition 3.3. We define the elliptic reconstruction G Hlifil) to be the unique 
solution of the elliptic problem 


B{t’‘;w\v) = {A\v) Vn G H^in). 


Remark 3.3. The dG discretisation of the equation in Definition 3^ is to find 
w\ G such that 


B(r-,wlvt) + K^(wlvt) = (A\vt) Vvi e Vj, 


Using the definition of A^, we obtain the equality 








Therefore = u\ and so B(t^; v) 


can be estimated using Theorem 


2.1 


At each time step fc, we decompose the dG solution u\ into a conforming part 
^0 (^) ^ ^ non-conforming part E such that c + d- 

Given t G we (re)de£ne Uh{t) to be the linear interpolant with respect 

to t of the values and viz., 


Uhit) := lk{t)ul + 


where {4,4+1} denotes the standard linear Lagrange interpolation basis dehned 
on the interval [t^,t^~^^~\. We dehne Uh^dt) and Uh^ifi) analogously. We can then 
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decompose the error e := u — Uh = ec — Uh,d where Cc ■= u — Uh,c- It will also be 
useful to dehne ■.= — u^. 

Remark 3.4. Aside from the numerical method itself and the bilinear forms 
appearing in Definition \3 . t^ then given t G bilinear forms appearing in 

the error analysis are assumed to take place over the union triangulation 
The norms ||| • ||| and \ ■ \a (and thus by extension || ■ ||*j are all evaluated over 
the union triangulation. 

Lemma 3.1. For almost any t G (t^,t^~^^] nie have 




Proof. This follows from (3.2). 


□ 


Before proving the a posteriori bounds for the fully-discrete method, we introduce 
the error estimators. We begin by dehning the spatial estimator, rjs, by 


n—1 


n—1 


Vl ■■= l|e(0)|P + 7^ + Vk,j+i) + 


j=0 


j=0 


'n—1 


n—1 


{ I + + 1 ^i+wl4,i+l ( > 

0=0 / j=0 


where 




^S2J + 1 


J ■= ~ 11"^' + - a' ■ ~ 11 I \Ihe) 

K&Qo E&SiC^) 

e&s{cA ^ ^ Ee£™‘(C-’') 

yi+i _ ji+i yt+i _|_ ^h~ 

^ Ej+l 


JT 

e 


hkj 


^ E 

A'eC-’uc^+i 


LAK) 


Vkj+i 


EesicA 

■= E - 

£;g£:(C-^'uc^+i) 


r 7+1 7 

r < -< 

Ej+l 


LAE) 
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The time (or temporal) estimator, r]T, is given by 


n—1 

/. ^Ti,j +1 dt+min 

j=0 dto 





where 


,i+i : = ^ ^ 11 ^i+i - a) + Ij (a^' - a) < 11 , 

A'eC'’’uC''+^ 


+ Ij+i - 6 - V ■ a^'+i + V ■ a) 


u 


,i+i I h 
h \\L'^{Ky 


Theorem 3.3. The error of the fully-discrete method (3.5) satisfies the bound 


< 


\IvI + Vt- 


Proof. From Lemma 3.1 and Dehnition 3.3 we have 

+ B v) + ^ , 


(3.22) 


which upon adding and subtracting (1^ ( 


, n) and using Remark 


3.2 


gives 


(1?’^) = {f - + lk{A’"^^ - A^),v) -B{t;uh,v) 

+ B{t'^+^-u'l+\v)+B{t^+^-e^+\v)-lkB{t^+^-w^+\v) (3.23) 

+ 45(t", v) + (^f+^ - ^ • 

Finally, we add and subtract lkB(t^~^^; and lkB(t^; u^,v) to obtain the 
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primary error equation, viz., 

= {f - + lk{A’"^^ - A’^),v) -B{t;uh,v) 


(3.24) 


+ hB{t^-,e\v) + (^f+^- 


4^+1 

dt 




By combining terms, using the definition of the bilinear form B and the 
Cauchy-Schwarz inequality; the first four terms give rise to the time estimator: 


(/ - + 4(41"+' - A'^),v) + 4+i5(t"+';n^+',n) + hB{t^-ut,v) 

- B[t;uh,v) < ?7Ti,fc+i||4||| + ?7T2,fe+i|4ll- 


(3.25) 


The final term can be rewritten using Remark 3.2 then bounded using Theorem 


2.4 and the Cauchy-Schwarz inequality, viz.. 




(3.26) 


< 


h52,fc+l|lTl 


For the remaining terms, we use the elliptic reconstruction property together with 
Theorem 12. II to conclude that 


4+iR(f''++e^+\n) + 4R(f'';e*',n) < (4+ih5i,fc+i + 4hSi,fc)||4l 


(3.27) 


Combining the above, setting v = and using (2.8), (2.9), the Cauchy-Schwarz 
inequality and Young’s inequality yields 


d 


') + 


< 


4+l%i,fc+l + 4hSi,fc + hS2,A:+l + 11 \'^hA 11 + \Uh,d\A 


(3.28) 


+ hS4,A:+l|4c|| + hri,fc+l + hT 2 ,fc+l||ec| 


The proof then follows from Theorem 3T and by employing a bounding strategy 
identical to that used in Theorem 13.21 □ 


Remark 3.5. The spatial estimator is expected to be asymptotically robust with 
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respect to e as the predominant terms are the same as in the elliptic case. For 
the pre-asymptotic case, one would need to work in stronger norms to achieve 
theoretical robustness. We note, however, that in all the numerical experiments 
below, the adaptive algorithm, implementing the estimators presented here, was 
able to arrive to guasi-optimal space-time mesh modifications. The temporal error 
estimator is also expected to be asymptotically robust with respect to e as the 
temporal data approximation error terms are all order two in time and the only 


order one temporal term is a difference of derivatives (from Remark 3.2) which is 
anticipated to be independent of e in the asymptotic regime. 

Remark 3.6. The use of elliptic reconstruction is not essential to the proof 


of Theorem 13.2\ and Theorem \3.3^ it is possible to derive the residual based a 
posteriori bounds directly albeit at the cost of a lengthier calculation. The advantage 
of using elliptic reconstruction in the proof lies in the fact that the space estimator 
can be easily modified to accommodate non-residual based elliptic error estimators. 
This, in turn, may offer improvements in robustness with respect to the Peclet 
number cf. m- 


3.4 An adaptive algorithm 

An adaptive algorithm is a computational procedure that seeks to use an error 
estimator, p, to try and minimise the error in some norm by appropriately reducing 
the discretisation parameters. For elliptic problems, such adaptive algorithms are 
relatively straightforward: a hnite element solution is calculated on an initial mesh 
and its estimator evaluated then the regions of the mesh where the estimator is 
largest are targeted for refinement by the adaptive algorithm and the hnite element 
solution is recalculated on this new mesh; the algorithm continues in this fashion 
until the error estimator is below a given tolerance. For parabolic problems, the 
design of adaptive algorithms is far more challenging because it is unclear how 
the spatial and temporal components of the error estimator should be utilised. 

The algorithms currently in the literature [301 001 0U 02l 03l 011 EHl EHl EQI 
|96] consist of an initial condition tolerance to control rehnement of the coarse 
input mesh, a spatial refinement tolerance to control mesh rehnement, a spatial 
coarsening tolerance to control mesh coarsening and a temporal tolerance to control 
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the length of each time interval. Typically, these algorithms focus on the use 
of these individual tolerances to force the error estimator below a given global 
tolerance. However, it is not necessarily clear that this is the correct choice. 
Indeed, proving that an adaptive algorithm will terminate with the total estimator 
below a tolerance is not the same as showing that it produces a quasi-optimal 
distribution of time steps and mesh parameters. 

We shall introduce a new adaptive algorithm, based on that given in [30] , with 
a different emphasis on the use of tolerances and we will show numerically that 
our adaptive algorithm reduces the error estimator that it utilises at the optimal 
rate with respect to the mesh parameters and the total number of time steps. 
The pseudocode for our algorithm is given in Algorithm 3.1 and is based on using 


different parts of the a posteriori estimator from Theorem 3.3 to drive space-time 
adaptivity. It is useful to provide heuristic justihcations for the approaches taken 
in our adaptive algorithm and to compare our adaptive algorithm to similar 
algorithms already in the literature: 


• As in |30|, our adaptive algorithm uses the dominant term in the space 
estimator, rjs^j+i, to control mesh rehnement. All elements on which 

is larger than the spatial refinement threshold stol+ are targeted for 
rehnement by the adaptive algorithm. 

• Most algorithms in the literature conduct mesh coarsening through the 

term T]s 2 ,j+i (or equivalent) which is often referred to as a mesh-change 
indicator — this is because such a term is non-zero only on elements that 
have been subject to coarsening. This approach, however, comes with a 
serious disadvantage — spatially one order higher than which 

means such algorithms tend to be too conservative with regards to mesh 
coarsening. By contrast, our adaptive algorithm uses the dominant term 
in the space estimator, 775 ^ j+i, to control mesh coarsening as well as mesh 
rehnement. In particular, all elements on which j+i is smaller than the 
spatial coarsening threshold stol“ are hagged for coarsening by the adaptive 
algorithm. 

• The nature of the time estimator, rjT, makes it inconvenient to use as a 
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Algorithm 3.1 Space-time adaptivity 
1: Input: e, a, b, f, Uq, T, 12, n, 7, ttol, stol’*', stol~. 

2: Set Ti, Tn = T/n. 

3: Calculate 
4: Calculate u\ from 

5: while 17 !’ 1 > ttol OR max 77 ^^ i\k > stol’*' do 

6 : Modify C,^ by refining all elements such that i\k > stol+ and 

coarsening all elements such that 77 !^ i\k < stol“. 

7: if 7)2 > ttol then 

8: n n + 1. 

9. Tjj Tyi—\i T 3 T 2 - 

10: T2 = Ti/2. 

11: T\ i — T\/2. 

12: end if 

13: Calculate u/. 

14: Calculate u\ from 

15: end while 

16: Set j = 1, time = ri. 

17: while time < T do 
18: Calculate from u\. 

19: while f/r j+i ^ ttol do 

20: if f/r j+i > ttol then 

21: n n + 1. 

22 : Tfj Tn—li •••1 hj+S Tj-\-2' 

23: Tj+2 = Tj+i/2. 

24: Tj+i ^ Tj^i/2. 

25: end if 

26: Calculate from 

27: end while 

28: Create from (t Py rehning all elements such that VsiJ+iIk > stol"*" 

and coarsening all elements such that < stol“. 

29: Calculate from 

30: time time -|- Tj+i. 

31: j^j + 1. 

32: end while 
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temporal refinement indicator so we define 57 tj+i given by 






Vtj+i 


7]Ti,j+idt + mm{aT,T} / VT2,j+idt, 


(3.29) 


the sum of which bounds rj'^. Our temporal strategy consists of continually 
halhng the time step length until fir j+i is smaller than the temporal threshold 
ttol which is in contrast to the approach taken in [20] where they instead 


seek to force the integrand in (3.29) below a given tolerance. We will show 
numerically that our approach yields a quasi-optimal distribution of time 
steps. 


• The hnal major difference between our proposed algorithm and those in 
the literature lies in how the coarse input mesh is dealt with. Typically, 
mesh rehnement is carried out on the coarse mesh until the initial condition 
estimator, | |e(0) 11, is smaller than a given tolerance. Such a term is, however, 
spatially one order higher than r]Si,i which means using it for conhguration of 
the input mesh leads to mass mesh rehnement during the hrst time step and 
such a large amount of mesh change during one time step can destabilise the 
numerical scheme. Therefore, our adaptive algorithm continues to rehne the 
coarse input mesh until ^ is smaller than the spatial refinement threshold 
stol+ on every element. 

Remark 3.7. Mesh modification must be done very carefully to ensure that the 
numerical solution does not degrade ITR Specifically, the spatial refinement 
threshold stoT*" needs to he chosen sufficiently small in comparison to the temporal 
threshold ttol. The spatial coarsening threshold stol“ also needs to be chosen 
sufficiently small in comparison to the spatial refinement threshold stol’*' in order 
to avoid unnecessary refine and coarsen loops. 

Remark 3.8. We stress that the algorithm does not necessarily produce a 
monotonically decreasing time step distribution from t) to T. Indeed, the algorithm 
starts with an initial eguispaced subdivision of [0,T] into n time intervals, which is 
then, possibly, locally bisected based on ttol. For instance, if the solution reaches 
a smoothly varying steady state, the algorithm will retain the original (coarse) 
time step length ofT/n during the final stages of the computation. 
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3.5 Numerical experiments 

We shall numerically investigate the presented a posteriori bounds and the 
performance of the adaptive algorithm through an implementation based on the 
deal. II hnite element library m- All the numerical experiments have been 
performed using the high performance computing facility ALICE at the University 
of Leicester. 

In order to discuss the numerical results, we need some additional dehnitions. 
We shall begin by extending our notion of the effectivity index. Let the maximum 
meshsize be given hy h := max max hx and the time largest step length be given 

0<k<n 

by At := max Tfc. We then dehne the spatial effectivity index by 

l<k<n 


spatial effectivity index := lim 


V 


At->o I |e| 

whereas the temporal effectivity index is given by 


temporal effectivity index := lim 




h^o e 


These notions give us a way of measuring the contribution of the spatial and 


temporal estimators to the constants in (2.11)-(2.12). In particular, we can 


assess whether specihc parts of the estimator are robust or not. The spatial 
effectivity indices can be observed in practise by choosing a very small temporal 
threshold while the temporal effectivity indices can be observed through use of a 
high polynomial degree and/or sufficiently hne spatial mesh. 

We also need a notion of the average number of degrees of freedom so we can 
discuss spatial convergence rates. If the total number of degrees of freedom on 
the union mesh Cf U is denoted by then the weighted degrees of freedom 
of the problem is given by 


n—1 


Weighted Average DoFs := — E 


3=0 


In all examples presented below, unless otherwise stated, we use polynomials 
of degree two and an initial 4x4 uniform quadrilateral mesh. We also set the 
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spatial coarsening parameter to stol“ = 0.001 ♦stol"'". Finally, unmarked lines in 
convergence plots represent the theoretically expected optimal rate of convergence 
for reference purposes. 


3.5.1 Example 1 

Let = (0,1)^, a = (1,1)"^, 6 = 0, Mq = 0, T = 10 and select the function / so 


that the exact solution to problem (3.2) is given by 

,(x-l)/£ _ 


u 


ix,y,t) = (1 - e *) 


3-l/e 


+ X — 1 




3-l/e 


+ y-l 


The solution exhibits boundary layers at the outflow boundary of the domain of 
width 0{e) as well as a temporal boundary layer. 

We begin by hxing a temporal threshold that produces enough time steps so 
that the temporal contribution to the error is very small in comparison to the 
spatial contribution. The spatial threshold is then gradually reduced to observe 


the spatial effectivity indices for this problem which are given in Figure 3.1 


Optimal rates of convergence are observed with respect to the weighted average 
degrees of freedom for both the estimator and the error but are omitted in this 
example. As shown, the effectivity indices are bounded asymptotically and remain 
between hve and ten for the different values of S] these are directly comparable 
to those observed in [97] for the stationary problem. 

In order to study the temporal effectivity indices for this problem any boundary 
layers must be fully resolved so that the spatial error is dominated by the temporal 
one. To this end, we use a high polynomial degree and specially constructed 
anisotropic meshes in order to ensure that the spatial error is sufficiently small. 
The temporal threshold is then reduced to observe the temporal effectivity indices 
of the problem which are given in Figure |3.1[ Optimal order is observed in both 
the estimator and the error and the effectivity indices remain bounded between 
four and seven for all values of £. 

The presence of a temporal boundary layer in the solution motivates a 
comparison between adaptive and uniform time-stepping. To this end, a sufficiently 
small spatial threshold is chosen so that the spatial contribution to the error is 
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Spatial Effectivity Indices Temporal Effectivity Indices 




Figure 3.1: Example 1: Spatial and temporal effectivity indices. 




Figure 3.2: Example 1: Temporal error comparison under adaptive and uniform 
time-stepping for e = 1 and e = 10“^. 


small and then the temporal threshold is decreased and the results are compared 


to just using uniform time-stepping. The results given in Figure |3.2| show that 
the temporal strategy of the adaptive algorithm minimises the temporal portion 
of the error better than just using uniform time stepping. 


3.5.2 Example 2 

We set n = (—1,1)^, a = (1,1)'^, 6 = 1,/ = sm{5t)xy, mq = 0 and T = 27r. The 
solution exhibits layers of width 0{e) in the proximity of the outflow boundary 
and is oscillatory in time. The sharpness of the boundary layers depend on time, 
thus making this a good test of the ability of the algorithm to add and remove 
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Rates of the Spatial Estimator 


Rates of the Temporal Estimator 




Weighted Average DoFs Time Steps 


Figure 3.3: Example 2: Spatial and temporal rates. 

DoFS vs Time fore = 1 ^^^4 DoFS vs Time for e = 10"^ 




Figure 3.4: Example 2: DoFS vs Time for e = 1 and e = 10 


degrees of freedom. 

As in Example 1, we begin by fixing a temporal threshold while decreasing the 
spatial threshold to observe the rates of convergence for the space estimator. We 
then set a spatial threshold small enough to resolve any boundary layers, while 
reducing the temporal threshold to observe the rates of the time estimator. The 


results are displayed in Figure [T3| Optimal rates of convergence are observed for 
both the space and time estimators. 

To assess the mesh change driven by the adaptive algorithm we also plot the 
individual degrees of freedom on each mesh against time for a given spatial and 
temporal threshold. The results are given in Figure 3.4 We observe that the 


adaptive algorithm is adding and removing degrees of freedom at a rate that is 
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Temporal Effectivity Indices 


Rate of the Time Estimator for e = 10 ^ 




Figure 3.5: Example 3: Temporal effectivity indices and the rate of the time 
estimator for £ = 10“^. 

in accordance with the oscillating nature of the solution driven by the sinusoidal 
forcing /. 


3.5.3 Example 3 


Let n = (-2, 2)2, T = 2%, a = {y, -x)^, 6 = 0, / = 0 and mq = 

The PDE convects the initial two dimensional Gaussian prohle along the circular 
wind while diffusing it at a rate depending upon e. In particular, provided the 


error at the boundary is sufficiently small, the exact solution to problem (3.2) is 
given by 


u{x,y,t) = 


1 + 256et 


exp 


64(a; — 0.5 cos(t))' 
1 + 256et 


exp 


64(2/ + 0.5sin(t))" 
1 + 256et 


To observe the temporal effectivity indices and temporal rates of the problem 
we hrst £x a spatial threshold so that the spatial contribution to the error is small 
and then reduce the temporal threshold; the results given in Figure [3^ show that 
the temporal effectivity indices are bounded and remain between one and eight 
for all values of £ and that the optimal rate of convergence is achieved by both 
the error and the estimator. Some meshes at various time steps produced by the 
algorithm for e = 10“^ are displayed in Figure 3.6 and show that the adaptive 


algorithm is adding and removing degrees of freedom efficiently. 
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Grid at t = 0 


Grid at t = Jt/2 



Grid at t = 7t 



Grid at t = 37i/2 



Figure 3.6: Example 3: Grid snapshots. 


3.5.4 Example 4 

Let G = (0,1)^, a = (sin(t),cos(t))^, 6 = 0, / = 1, mq = 0 and T = 27r. The 
nature of the solution is rather uniform in time but spatially the solution possesses 
a moving boundary layer of width 0{e) that is driven by the changing nature of 
the inflow and outflow boundaries. Therefore, this example is well suited to testing 
the ability of the algorithm to adapt the grid to this moving boundary layer. Grids 


at various times are shown in Figure 3.8 for £ = 10 


^-2 


As in previous examples, we £x a small temporal threshold and then reduce 
the spatial threshold to observe the rates of the space estimator. Again, we also £x 
a spatial threshold small enough to ensure that all boundary layers are sufficiently 
resolved and then reduce the temporal threshold to observe the rates of the time 


estimator. These results are given in Figure 3.7 


Optimal spatial and temporal rates of convergence are observed and the grids 
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Space Estimator 




Figure 3.7: Example 4: Spatial and temporal rates. 




Grid at t = 3it/2 




Figure 3.8: Example 4: Grid snapshots. 
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produced for e: = 10“^ clearly show that the adaptive algorithm is picking up the 
boundary layers as they move around the domain and that unneeded degrees of 
freedom are not retained. 

3.6 Conclusions 

An a posteriori error estimator for the discontinuous Galerkin spatial discretisation 
of a non-st at ionary linear convection-diffusion equation was derived. The numerical 
examples presented clearly indicate that the error estimator is practical and the 
respective space-time adaptive algorithm works well for the studied problems. As 
predicted, the spatial effecitivity indices are in an identical range to those observed 
in Ea and the spatial part of the error estimator appears to be asymptotically 
robust with respect to e. Furthermore, the temporal effectivity indices of the 
studied problems are substantially smaller than those seen in [SS] for the heat 
problem and may even be fully robust with respect to the norm || ■ ||*. 
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Chapter 4 

A posteriori error estimation 
and blow-up detection for 
nonlinear ODEs 


4.1 Blow-up in nonlinear ODEs 


This section is devoted to the numerical approximation of the ODE 

du 

Tt = 

m( 0) = Mo, 


(4.1) 


where / is a Lipschitz continuous function and uq > 0. We say that (4.1) exhibits 
blow-up if the solution u has the property that 


limsup \u{t)\ = oo, 

t^T* 

for some T* > 0. The value T* is referred to as the blow-up time and if T* < oo 
we say the ODE exhibits finite time blow-up. Throughout the rest of this chapter, 


it is assumed that (4.1) exhibits hnite time blow-up. Regarding the numerical 


approximation of (4.1), the following questions are of interest: 


Using a simple time stepping scheme to approximate (4.1), can we construct 


a residual based a posteriori error estimator for the given method? 
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• If the answer to the above is yes, is it possible to use the estimator to drive 
an adaptive algorithm that successfully converges to T*? 

• If the adaptive algorithm does converge to T*, is it possible to numerically 
quantify how quickly different time stepping schemes converge to T*1 


In attempting to answer the above questions, there are some papers in the 


literature of particular interest. In |98], the authors approximate (4.1) using a 
6'-method along with a temporal rescaling of the ODE; this modihes a distribution 
of uniform time steps so that they better match the blow-up behaviour of the 
numerical solution. Their numerical solution displays the same asymptotic 
behaviour as the exact solution. 

Similarly to [98], in (GOj the authors also transform the ODE but through an 
arc length transformation. They use a forward Euler method to approximate the 
transformed equation and they show that their adaptive algorithm, which is based 
on their transformation plus a tolerance controlled ODE integrator, converges 
towards the blow-up time linearly with respect to the total number of time steps. 
However, in contrast to [98], they restrict themselves to the case f{u) = u^, p > 1. 

Finally, in [US] they prove existence results for numerical approximations to 


(4.1) provided that the time step lengths are sufficiently small and that the 
nonlinearity satishes a polynomial growth condition. For the specihc case 
f{u) = u^, p > 1, they show that a certain selection of time step lengths yields 
approach to the blow-up time. 


4.2 An a posteriori error estimator 


For simplicity, we restrict our attention to 


/(“) = 

3=0 


(4.2) 


where p > 2 is some positive integer and the coefficients satisfy Cj > 0 (cp > 0) 


so that the problem is guarenteed to blow-up. In order to approximate (4.1), we 


shall use a generic one-step scheme with right-hand side fh that approximates /. 
That is, we set = Uq and for /c > 0 with some time step length r^+i, we search 


50 





for such that 


_ „,k 

“ft “ft 


'Tk+l 




(4.3) 


We also recursively dehne our time \= + Tk+i with := 0. In order to 

discuss the error of different time stepping schemes, we need to describe Uh on the 
interior of the intervals . Thus, given t G (t^, , we dehne Uhit) to be 

the linear interpolant with respect to t of the values and viz.. 


Uh{t) := lk{t)ul + 


It is now possible to construct an error equation for (4.3) by subtracting (4.3) 


from (4.1). Then, dehning e := u — Uh, we obtain 




(4.4) 


Adding and subtracting f{uh) to the right of (4.4) and dehning the residual 


r]k+i := f{uh) — /ft(M^, we obtain the error equation 


* = + /'(«,.)e + ± 


dt 


i=2 


3'- 


(4.5) 


where denotes the order j partial derivative of / with respect to u. In order to 


derive a usable error estimator from (4.5), we make use of Gronwall’s inequality. 


Application of Gronwall’s inequality to (4.5) for t G yields 


|e(t)| < Hk+i{t)Gk+i4>k+i, 


(4.6) 


where 


p pt 

Hk+i{t) := exp ( ^ 


Ifk 


f^^Kuh) 


3'- 


P ^ ds I , 


Gk+i 


i=2 

exp 1 / \f{uh)\ds I , 

Itk 


0fc+i := |e(f^) I + / \r]k+i\ds. 


l^k 
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Note that this is not truly a posteriori yet due to the presence of Hk+i- order 
to amend this, we use a local continuation argument in the spirit of [131 EH El]- 
To this end, dehne the set 


4+1 := t e 


max |e(s)| < 4+iGfc+i0fc+i 4 


where 4+i > 1 is a parameter to be chosen. Obviously G 4+i so 4+i is 
non-empty and bounded. We denote the maximal value of t that belongs to 4+i 


by t* and we assume that t* < From (4.6), we have that 


max \e{s)\ < Hk+i{t*)Gk+i(j)k+i- 


(4.7) 


By the dehnition of the set 4+1, we have 


0=2 

p 

0=2 


e(5)|)' ' 1 


fU 



J J 

/ 

Jtk 

f^^Kuh) 

j! 


ds 


ds . 


(4.8) 


Therefore, 


max |e(s)| < Gk+iA+iexp | ^i+i^klWk+\ / 


f^^Kuh) 


J'- 


ds . (4.9) 


Now, suppose that the upper bound in (4.9) is bounded strictly from above by 
the upper bound of the set 4+i, viz.. 


,tk + l 


Gfc+i0fe+i exp 

0=2 


k+l^k+l 


Ifk 




ds j < dk+iGk+i4>k+i, (4.10) 


or equivalently. 


exp ( J2^i+\Gi+\(l>k+\ 
0=2 




<tk 


f^^Kuh) 


J'- 


ds j <c (5/c-(-i7 


(4.11) 
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then t* cannot be the maximal value of t that belongs to Ik+i because we just 
showed max |e(s)| satishes a bound strictly less than that assumed in the set 

se[4'=,t*] 

Jfc+i — a contradiction. Therefore, provided (4.11) is satished, Ik+i = 
and we have our desired error bound once we select 5fc+i. Given that we wish to 
construct the best bound possible, we seek to minimise (4.11). Taking the limit 
we can, in fact, just select Sk+i to be the minimiser of 


1=2 


k+1 



f^^Kuh) 


j'- 


ds - log(4+i) = 0, 5k+i > 1. (4.12) 


Therefore, providing the solution to (4.12) exists, we have the following error 
bound 

|e(t)| < 4+iGfc+i0fc+i- (4.13) 


|e(f^+GI< max 


Remark 4.1. The term (j)k+i can be redefined with |e(t^)| estimated using the 
error estimator from the previous time step without any loss of generality to the 
argument presented giving us a recursive procedure for estimating the error. 


Remark 4.2. In practice, the solution to (4.12) is approximated using a Newton 
method. 


A natural question that arises is whether or not (4.12) can be satished practically 
close to the blow-up time. With the aid of the next lemma, we state a precise 


condition on the time step lengths t^+i which indeed ensures that (4.12) has a 
root > 1. 

p p 

Lemma 4.1. If''^jCje^ < 1 then s{x) = ''^CjX^ — log(a;) with Cj > 0, j = 1, 


1=1 


1=1 


..., p, p G N has a root in (1, -|-oo). 
Proof. See Lemma 2.2 in [24] . 


□ 


The above lemma gives a sufficient condition on when (4.12) can be satished. 


In particular, condition (4.12) can always be made to be satished provided that 


the time step length Tk+i is chosen such that 

J- Jtk 


\f^^\uh)\ds < 1 . 


1=2 


53 














It is useful to discuss in heuristic terms what 6k+i, Gk+i and (j)k+i in (4.13) 


represent. Obviously (pk+i is an approximation to the error on each time interval 
but what about Gk+i and 5^+1? Clearly both Gk+i and 6k+i are 
accumulation factors that represent the contribution of blow-up to the error 
estimator in some way. To gain some insight on these multiplicative terms, 
consider f{u) = u^, p > 1. Through separation of variables, the solution to 


(4.1) is given by 


u[ 


1 

1-p 


Now, suppose that Uh 


i(t) = {ul P + {l-p)t) 
u and consider Gk+i{u) given by 


Gk+i{u) = exp I / f\u) ds I = exp ( / pvF~^ ds ) . 


Substitution of the exact solution yields 


= exp 



Ga:+i(m) = exp I / p{ul ^ + {I - p)t) ^ ds 


ul ^ - 1 - (1 - 

ul~^ + (1 - p)t^ 


u 


i-p + (l_p)tfc+i\ 


1-p 


Wo ^ + (1 - p)t’" ) 


So for f{u) = uP, Gk+i{u) measures the blow-up rate of the exact solution on 
the interval and thus 6k+iGk+i can be viewed as the blow-up rate of the 

numerical solution. As we performed a Taylor expansion in our error analysis, we 
infer that Gk+i is the linearised numerical blow-up rate on the interval 
and Sk+i is the higher order part of the numerical blow-up rate on the interval 


. With these notions, another way of viewing (4.11) is that the numerical 
solution ceases to be valid once the (approximate) higher order terms from the 
Taylor expansion start to become dominant in the error estimator. 
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4.3 Adaptivity and convergence towards the blow-up time 


Using our knowledge of algorithms utilising a posteriori error estimators for linear 
problems, we propose Algorithm 4.1 for advancing towards the blow-up time. The 
basic idea behind the algorithm is to half the time step length and recompute the 
solution until the residual is below a given input threshold tol. The algorithm 
then advances by using the previous (now hxed) time step length as a reference 
to compute the next approximation. The algorithm continues in this way until 


(4.12) no longer has a solution; the algorithm then terminates and outputs the 
total number of time steps N and the hnal time T. 


Algorithm 4.1 ODE Algorithm 1 
1: Input: /, fh, Uq, Ti, tol. 

2: Calculate u\ from u^. 

3: while / |? 7 i| ds > tol do 

JtO 

4: Ti ^ ri/2. 

5: Calculate u\ from 

6: end while 
7: Calculate (5i. 

8: Set k = 0. 

9: while Sk+i exists do 
10: k k + 1. 

11- Ifc+l 'T~k- 

12: Calculate from 

13: while / |? 7 fc+i| ds > tol do 

14: Tfc+i t— rfc+i/2. 

15: Calculate from u\. 

16: end while 

17: Calculate d^+i. 

18: end while 
19: Output: k, 


Assuming that the adaptive algorithm outputs successfully, we wish to observe 
the order with which the adaptive algorithm approaches the blow-up time. To 
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this end, we define the A function 


A(tol,iV) := |T* -T(tol,Ar)| , 


where T* is the blow-up time of problem (4.1). It is conjectured that 


X{tol,N) oc iV-", 

where r is the order with which the adaptive algorithm approaches the blow-up 
time. An educated guess would be that r is the same as the order of the method 
that we choose to use. 

In order to gain some insight on how A converges, we apply Algorithm 4.1 to 


problem (4.1) with f{u) = for p = 2, 3 and m( 0) = 1 under the following time 


stepping schemes 


Explicit Euler 


Implicit Euler 


Improved Euler 



) = /K)^ 

)=/(«n> 


) = + /K + T-fc+i/K)))- 


The approximate rates of A under Algorithm 4.1 are given in Table 4.1 


Table 4.1: ODE 

Algorithm 

1 Results 

Method 

p = 2 

p = 3 

Implicit Euler 

r 0.66 

r 0.79 

Explicit Euler 

r ^ 1.35 

r ^ 1.60 

Improved Euler 

r ^ 1.2 

r 1.48 


The first question that arises is why is the explicit Euler method significantly 
better than the implicit Euler method? The answer lies in the way in which 
we have derived the error estimator. The numerical solution from the explicit 
Euler method always underestimates the true solution u ^8]; this means 5k+i is 
correcting for the fact that Gk+i is underestimating the true blow-up rate — our 
error bound is very tight and this explains the high convergence rate of A. For 
the implicit Euler method, Gk+i overestimates the true blow-up rate [98] meaning 
we obtain nothing “extra” from the error analysis in the way that we do for the 
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explicit Euler method. 

The second question is why is improved Euler worse than explicit Euler? 
Indeed, one would expect a faster approach to the blow-up time with a higher 
order method. The reason for this lies in a fault with the proposed adaptive 
algorithm. Indeed, the threshold approach taken to reducing our time step length 
is good for linear problems. However, we have neglected the presence of Gk+i in 
our error estimator; this factor tells us that each successive interval matters less 
to the error estimator than previous intervals meaning we need to increase tol 
on each interval. Thus, we propose Algorithm 4.2. 


Algorithm 4.2 ODE Algorithm 2 
1: Input: /, fh, uo, Ti, tol. 

2: Calculate u\ from 

3: while / It^iI ds > tol do 

JtO 

4: Ti ■(— Ti/2. 

5: Calculate u\ from u^. 

6: end while 
7: Calculate (5i. 

8: tol = Gi * tol. 

9: Set k = 0. 

10: while 6k+i exists do 
11: fc <(— fc -|- 1. 

12: T~k+1 'T~k- 

13: Calculate from 

14: while / |? 7 fc_|_i| ds > tol do 

Jtk 

15: Tk+i ■<— Tfc_|_i/2. 

16: Calculate from 

17: end while 

18: Calculate 6k+i. 

19: tol = Gk+i * tol. 

20: end while 
21: Output: k, 


The rates of A under Algorithm 4.2 are given in Table 4^ and a more detailed 
comparison of the convergence of A for the different algorithms under the various 


time stepping schemes is given in Figures 4.1, 4.2 and 4.3 These results show 
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Figure 4.1: Rates of A for explicit Euler under ODE Algorithms 1 and 2. 




Figure 4.2: Rates of A for implicit Euler under ODE Algorithms 1 and 2. 




Figure 4.3: Rates of A for improved Euler under ODE Algorithms 1 and 2. 
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that under Algorithm 4.2 we have recovered our conjectured rates for the function 
A with a slight bonus rate for the explicit Euler method. Note that for p = 3, 
Algorithm 4.1 converges faster than Algorithm 4.2 for the explicit Euler method; 
it is unknown why this is. 


Table 4.2: ODE Algorithm 2 Results 


Method 

p = 2 

p = 3 

Implicit Euler 

r ^ 1.00 

r 1.00 

Explicit Euler 

r ^ 1.45 

r 1.43 

Improved Euler 

r ^ 2.03 

r ^ 2.03 


4.4 Conclusions 


We derived an a posteriori error estimator for a class of nonlinear ODEs exhibiting 
blow-up and applied the estimator in two different adaptive algorithms to try and 
approximate the blow-up time. Both algorithms converged to the blow-up time in 
all test cases with Algorithm 4.2 outperforming Algorithm 4.1 in almost all test 
cases. In particular, we infer that explicit treatment of nonlinearities of the form 


(4.2) appears to be advantageous in the context of adaptive algorithms based on 


rigorous a posteriori bounds. 
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Chapter 5 

Adaptivity and blow-up detection 
for nonlinear non-stationary 
convection-diffusion problems 


5.1 Blow-up in semilinear PDEs 


For (non-fixed) T > 0, we consider the model problem of finding m : x (0, T] —)■ M 
such that 


— — eAu + a ■ Vm + f{u) = 0 

u = 0 

m (-, 0 ) = uo 


in X (0,T], 

on dil X (0, T], 
in 


(5.1) 


It is assumed that the reaction term f{u) is of the form f{u) = fo — v? although 
more general nonlinearities can be considered as discussed later in this chapter. 


We say that (5.1) exhibits blow-up if the solution u has the property that 


limsup ||M(t)||z,oo(j^) = cx), 

t^T* 


for some T* > 0. The value T* is referred to as the blow-up time and if T* < oo 


we say the PDF exhibits finite time blow-up. If (5.1) exhibits hnite time blow-up 


then we can describe the asymptotic spatial behaviour of the solution u through 


60 





Figure 5.1: Numerical approximation showing a typical solution profile near the 
blow-up time. 

the blow-up set, B, given by 


B := {x eQ \ 3{Xn, tn} C fix (0, T*), -)■ X, u{Xn, tn) Oo}. 

Elements of the blow-up set are referred to as blow-up points. The asymptotic 


spatial behaviour of the solution to (5.1) can be described through the blow-up 


set to be in one of two separate categories [ini|52]: 

• Point blow-up — B consists of a hnite number of points. 

• Regional blow-up — The one-dimensional Hausdorff measure of B is positive. 

For single point blow-up, the solution looks like a nascent delta function close 
to the blow-up time (see Figure multi point blow-up or regional blow-up can 
cause even more complicated spatial behaviour near the blow-up time. These 
demanding and complex spatial and temporal features make the numerical 


approximation of (5.1) and related PDEs very difficult and necessitates the 


development of adaptive hnite element methods. 
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ForT < T*, the weak form of Q reads: findw G L^{0,T-, H^{Q))nH^{0,T-, L\Q)) 
such that for almost every t G (0,T] we have 

( —,uj + + (/(t;M),u) = 0 Vu G Ffo(fi), (5.2) 

and the following assumptions are made on the PDF coefficients: uq G ipQ(f2), 

0 < £ < 1, a G [C(0, T; and /o G ^(O, T; L^(ff)). For simplicity, we 

assume that V ■ a = 0 and throughout the rest of this chapter, it is assumed that 


(5.2) exhibits hnite time blow-up. 


5.2 Space-time discretisation 

The numerical approximation of blow-up in nonlinear problems has been discussed 
in the literature. Solution prohles close to the blow-up time can be obtained 
through the rescaling algorithm of Berger and Kohn [T6| [85] or the MMPDE 
method [201 E5] • There is also work looking at the numerical approximation of 
blow-up in the nonlinear Schrodinger equation and its generalisations BlElllsol 
m [M]. Other numerical methods for approximating blow-up in a variety of 
different nonlinear PDFs can be found in |H1 1331 1331 SHI El]- Finally, in jlU2j . 
the author gives conditions that a numerical method must satisfy in order to 
asymptotically converge to the blow-up time. Most of these numerical methods 
rely on some form of theoretically justihed rescaling; however, there is no general 
theory to know whether the resulting numerical approximation is reasonable or 
not. In this chapter, based on the work contained in [23], we shall use a simple 


numerical scheme to approximate (5.2) and we will seek to perform our rescaling 


through rigorous a posteriori error estimates. 


We consider a full discretisation of problem (5.2) by using a hnite difference 


method to approximate the time derivative, taking the convection-diffusion terms 
implicitly for stability purposes and the nonlinear reaction term explicitly in 
view of the conclusions drawn in the previous chapter. To this end, consider 

n 

a subdivision of [0, T] into time intervals of lengths ri, ..., such that E = T 
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for some n > 1 then set := 0 and := Denote an initial triangulation 

j=i 

by and we further associate to each time step fc > 0 a triangulation which 
is assumed to have been obtained from by locally rehning and coarsening 
To each mesh we assign the hnite element space V/^ := Vh(C^) given 
by ( TSi and we also set := a(-,t^) and := for brevity. Finally, 

for t E (f^, , we let r denote the union of all edges in the mesh U 

The IMEX dG method then reads as follows. Set to be a projection of Uq 
onto 14 °. For k = 0, ..., n — 1, End G such that 


u 


fc+i 


u 




T'k+l 


for all G We shall take to be the orthogonal projection of 

onto 14 °, although other projections onto can also be used. 


Uo 


5.3 An a posteriori bound for the IMEX dG method 

The a posteriori error estimation of nonlinear parabolic problems has recently 
attracted attention for a variety of different PDEs [131 EH EH EH [32l [571 EH 11031 
110611107] . With regards to blow-up, a key result is by Kyza and Makridakis [75117^ 
wherein they produce an estimator for the error in the L°°{L°°) norm for a time 
semi-discrete approximation to the heat equation with polynomial nonlinearity. 
We shall use the ideas from [TSl ES] to derive an error estimator for the IMEX dG 


scheme (5.3) in the norm. A key novelty of the proof will be the use of a 

continuation argument for energy estimates, rather than the semigroup approach 
used in [THl [79] . 

Before we begin to construct our error bound, we require some additional 
notation. At each time step /c, we decompose the dG solution u\ into a conforming 
part u^c ^ n Vjf and a non-conforming part G Vjf such that = 

u’l c+u^ d- Further, given t G , we define Uh (t) to be the linear interpolant 

with respect to t of the values and viz.. 


Uh{t) := lk{t)u’^ + lk+i{t)u 


k-\-l 
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and we define Uh,c(t) and Uh,d(t) analogously. We can then decompose the error 
e := u — Uh = Cc — Uh,d where ec '■= u — Uh,c- If 'will also be useful to dehne the 


Lemma 5.1. Given t G then for any v G HQiil): 


elliptic error := 


— where w 


k ; 


is given as in Dehnition 3.3 


de 

~dG 



+ Bit] e, v) + if it] u) - fit] Uh),v) 


fif]Uh) 



Bit]Uh,v). 


Proof. This follows from (5.2) 


From Lemma [5. II we have 


□ 


+B{t]e,v) + {fit]u) - fit]Uh),v) = + f + 

+ -B{t]Uh,v) + {f - fit]Uh),v), 

which upon straightforward manipulation gives 


+B{t]e,v) + {fit]u) - fit]Uh),v) = + f + 

+ h+iB{t’^^^]e’^^\v) +hB{t^]e’^,v) - B{t]Uh,v) +h+iB{t^+^]ul+\v) 
+ hB (^^ ulv) + {f - f {t] Uh) + h {A^^^ - A ^), v ). 


We are now ready to state our a posteriori estimator. The hrst part of our 
estimator is the initial condition estimator, rji, given by 

hi ■= fl|e(0)|r+ ^i?||K]|lL2(ij)) • 

V E&£(C°) J 

Due to the nature of the error bound to be presented, it is easier to separate the 
remainder of the estimator into two parts. As in Chapter 3, a subscript S denotes 
parts of the estimator related to estimating space while a subscript T denotes 
parts of the estimator related to estimating time. In this way, for t G 
rjA is given by 


VA ■— lkVSi,k + lk+lVSi,k+l + VS2,k+l + VTi,k+l, 


64 









where 


VSi,k ■ = 


h \ \ , u . u u _ fc||2 


^ -^||A^ + 5AM^-a^-Vu 


I I r k fcl I |2 




^1 \l'^{k) ' i _/ g. 

£;e£(C'“) 




\L^(E) 


1/2 


E £ltK]"^ 

EeeiC^) 


L^E) 


+ J2 


fcl l|2 


E^£i^i{Qk) 


TlS2,k+l '■ — 


E 


hi 


K 


jk _ jk-\-l jk _ 


U 




h “ft 




'Tk+l 


^ii'eC'‘uc'“+^ 

?7Ti,fc+i := - a)M^+^ + 4(a*^ - a) 


lL2(E) 


L2(X) 


1/2 


U 


ftp 


while 1 /b is given by 


VB '■— VS3,k+l + h54,fc+l + VT2,k+l, 


where 


with 


1/2 


hS3,fc+i := 1 Y1 5^ ^I^EllMIlihE) 
^xec'=uc''+i Eci^is 

VSi,k+i ■= ( - - - 

\Ecr L ^fc+i J l^e)^ 

VT2,k+i ■= \ \f^ - /(4wft) + lk{A^^'^ - 


O-R ■- ‘^\\Uh\\L^{K) + IIMIlL°°(h-B)- 


Going back to (5.5), the hrst term on the right can be bonnded using Theorem 


2.4 and the Cauchy-Schwarz inequality, viz.. 




duk 

dt 


,v-irv 


(5.6) 


< VS2,k+i\\\v\ 
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The next two terms give rise to parts of the space estimator via Theorem 2.1 




(5.7) 


Using the dehnition of the bilinear form B and the Cauchy-Schwarz inequality, 
the hnal four terms give rise to the time estimator: 


- B{t]Uh,v) < r]T,,k+i\\\v\\\, 
+lk{A^^^ - A^),v) < VT2,k+i\\v\\. 


(5.8) 


Setting n = Cc in (5.5), using the results above along with coercivity of the bilinear 


form B and the Cauchy-Schwarz inequality we obtain 


+ llledir + - f{t-,Uh),ec) 


< 

rsj 


du 


h,d 


dt 


+ VT2,k+l I I Fc 


(5.9) 


+ {hVSi,k + h+lVSi,k+l + VS 2 ,k+l + hTi,A:+l)llkc||| + B(t; Uh,d-, Cc). 


Using continuity of the bilinear form B and Theorem 3.1 yields 


^^lledd + ll|ec||d + (/(^;«) - fit]Uh),ec) < (7754,^+1 + hT2,fc+i) Ike 


(5.10) 


+ ihVSi,k + lk+lVSi,k+l + VS2,k+l + hri,A:+l)|||ed||- 
Using Young’s inequality and the dehnition of our estimators, we conclude that 

^^Ikdd + ^llkdld + (/(Um) - f{t;Uh),ec) < ^hi + hslledl- (5-11) 
We must now deal with the nonlinear term. We begin by noting that 

(/(t; u) - /(t; Uh), ej = (/(t; + Uh) - /(t; Uh), ej = Ti + T 2 , (5.12) 


where 


^1 • {p‘'^h'^h,di {A'h,di^c)i 

T 2 ■ (2'U/iec, Cc) T (p‘(ic'^h,dy (®c)^c)' 

We can write the contributions to Ti elementwise and then use the Cauchy-Schwarz 
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inequality and Theorem 2.3 to conclude that 


1/2 


|r.| < ^ (2\\u^\\l °°{K) + \\Uh,d\\L'^[K)) ||ec 


< 


^A'eC''uc''+^ 

VS3,k+l I |Cc| I 


(5.13) 


To bound T2, we use Holder’s inequality along with Theorem 2.3 to conclude that 


1^2! ( 2 | |l°°(o) + l|['Wfe]||L°°(r)) Ike IP + Ik, 


c\\L3{n)- 


(5.14) 


Combining (5.11), (5.12), (5.13) and (5.14) we obtain 


with 


^IkclP + IlkcliP ^ + 2i/B|kc|| + 2crn|kc|P + 2|kc||i3(Q), 


■— 2| |m/i| kooj-Q) + C'||[M/i]||x,oo(r), 


(5.15) 


where C is a generic constant. We now note the Gagliardo-Nirenberg inequality 
which coupled with Young’s inequality yields 


1 

lkc|lL3(f2) Y iC||ec|| llVedl Y “||kc||| H Ike 

Combining these results with the generic constant C yields 


(5.16) 


^IkelP < Ct]\ + 2C'7/s|kc|| + 2crn|kc|p + 


2^-l| 


(5.17) 


We now need a way to deal with the norms appearing on the right-hand side. 


To do this, we use a variant of Gronwall’s inequality (see Theorem 2.13). In order 


to apply this to (5.17), we need to introduce some new notation. We dehne 


Gk+i := exp 


(^nds , 


Itk 

Hk+i{t) := exp ( [ \ \er.\\'^ ds 
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Then application of Theorem 2.13 to (5.17) for t G yields 




(5.18) 


where 


ife+i \ 1/2 ,fc+i 

0fc+i := ( \ \ec{t^)\\‘^ + C ri\ds\ +C r]Bds. 


We now need to remove the term iffc+i from (5.18) in order to construct a usable 


estimator. In order to do this, we use a continuation argument. To that end, we 
dehne the set 


h+i ■— {i e I ||ec||L°°(t'=,i;L2(Q)) < 6k+iGk+i(j)k+i}, 

where 5^+1 > 1 should be chosen as small as possible. We know that Ik+i is 
non-empty since £ Ik+l and obviously bounded. Denote the maximal value of 


t in Jfc+i by t* and assume that t* < then from (5.18) we have 


|ec||L°°(ffc,t*;L 2 (o)) < H(t*)Gk+i(l>k+i 

^ exp £ Tk+\\\£c\\ijOo(j^kk+l(t^k+l 

< exp (^K‘^£-^Tk+i5l^iG fc+i^fc+ijCfc+i^fc+i- 


(5.19) 


Now, suppose that 


exp {^K^e ^Tk+i^t+iGlkj^ift^^j^^Gk+ift^k+i < ^fc+iCfc+i^fc+i, (5.20) 


or equivalently. 


K £ < log(4+i), 


(5.21) 


then t* cannot be the maximal value of t in Ik+i because we just showed that 
||ec||L°°qfc,t*;L 2 (f 7 )) satishes a bound strictly less than that assumed in the set Ik+i 


— a contradiction. Therefore, providing (5.21) is satished, Ik+i = and 

we have our desired error bound once we select Sk+i- Taking the limit, we can 







select ^fc+i to be the minimiser of 


K‘^e - log(4+i) = 0, 4+1 > 1. (5.22) 


In order to obtain our error estimator, all that remains is to estimate (/>!. Application 


of Theorem 2.3 and the triangle inequality yields 


leMW" < MOW + \M0)\\^ < Cvl 


Therefore, if we (re)de£ne 0i to be 


\ 1/2 


')i ^ CT]j + C I rj\ds\ +C I T]Bds, 


I to 


I to 


then we have 

||£'c(t )|| 4 I l^-cl |L°°(t0,tl;L2(O)) 

where V’l := diGitpi. In the same way, if we (re)define 

\ 1/2 


„tk + l 
'.2 , / „2 


0fe+i ■=\'4’k + C Va 

Jt>^ ) 

i/’fc+i := 4+iGfc+i0fc+i, 


ds I +C 


tk+i 


Itk 


Vb ds, 


then we have 

||6c(t )|| 4 I l^'cl — '4^k+l- 

Hence, the following result holds. 


(5.23) 


(5.24) 


(5.25) 


Theorem 5.1. The error of the IMEX dG discretisation of problem (5.2) satisfies 

. 1/2 

|e||L-(0,T;L2(fl)) < V'n TeSSSUp 


0 <t<T 


\Ecr 


providing that the solution to (5.22) exists for all time steps. 


Proof. Follows from the above derivations, the triangle inequality and the bounds 
in Theorem 12.31 □ 
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The estimator produced above is suboptimal with respect to the mesh-size 
as it is only spatially optimal in the norm. It is possible to conduct 

a continuation argument for the norm rather than the norm if 

one desires a spatially optimal error estimator; this is stated for completeness 
in the theorem below. However, the resulting S equation was observed to be 
more restrictive with regards to how quickly the blow-up time is approached. For 


this reason, we opt to use the a posteriori error estimator of Theorem 5T in the 
adaptive algorithm introduced in the next section. 


Theorem 5.2. The error of the IMEX dG discretisation of problem (5.2) satisfies 

1/2 


f>T \ 1 /^ ^ 

E i/’fc -t- ess sup 


|e(r)|/+ / £||Ve|/* 


'0 


k=l 


0<t<T 


\Ecr 


Furthermore, close to the blow-up time where ||e(T)|| = | |e| |i,oo(o,T;L 2 (o)) we have 


\ 1/2 


< 


^ V'fc + ess sup [^hE\\[uh]\\l 2 ^^E)j . 


0<t<T 


k=l \Ecr 

where fjk, k = 1, n, is defined recursively with fjQ = Crji and 

( fk \ 1/2 

'ipl-i + C [ T]\ds + C [ plds] , 

Gk := exp(rfc/2) exp an , 

• ^kGk4^k'; 

provided that > 1 which is the smallest root of the equation 


Ke ^^‘^rj‘^5kGk(f)k - log(4) = 0, 


exists for all time steps. 


Proof. The proof is completely analogous to that of Theorem |5.1| and follows from 
(5.15) by conducting a continuation argument for the norm. □ 
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Remark 5.1. Although we considered a very simple nonlinearity, the continuation 
argument in this section can he modified to include any nonlinearity of the form 
f{u) = /o + fiU + f 2 U^ + fsu^- With a nonlinearity of this form, we would have to 
deal with the term | jed in the error eguation. From the Gagliardo-Nirenberg 
ineguality, we have 

WecWl^^n) < l|ec|n|Vedd. 

After application of GronwaWs ineguality, we have a term of the form 

( ^ 

exp I J 11 Ved d ds j , 

which can he hounded through a continuation argument that uses the norm. 

Any higher order nonlinearities could not he dealt with in this way. 

Remark 5.2. It is also worth noting that although we are primarily interested 
in looking at hlow-up problems, the estimators developed in this section are still 
perfectly valid if blow-up does not occur. 


5.4 An adaptive algorithm 


The a posteriori bounds presented in the previous section will be used to drive a 
space-time adaptive algorithm that is designed to approximate the blow-up time 


of problem (5.2). The pseudocode of this algorithm is given in Algorithm 5.1. 


As in Algorithm 3.1, both mesh rehnement and coarsening are driven by 
the term ? 75 i,fc+i. The size of the elemental contributions to ? 75 i,fc+i determines 
whether the elements are to be rehned, coarsened or neither depending on two 
spatial thresholds stol+ and stol". Similarly, rjT 2 ,k+i is used to drive temporal 
rehnement and coarsening subject to two temporal thresholds ttol+ and ttol” 
on each time interval. As in Algorithm 4.2, all spatial and temporal thresholds are 
increased by the factor Gk+i on the interval after the solution has been 

calculated. The algorithm then advances by using the previous (now hxed) time 
step length as a reference to compute the next approximation. The algorithm 


continues in this way until (5.22) no longer has a solution and the algorithm then 


terminates and outputs the total number of time steps, the hnal time and the 
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Algorithm 5.1 Space-time adaptivity 
1 : Input: e, a, /o, mq, 7 , ttol+, ttol^, stol'*', stol~. 

2 : Calculate 

3: Calculate u\ from u\. 

4: while / ? 7 ^ ^ ds > ttol^ OR max 77 ! i\k > stol"*" do 
JtO ’ ^ 

5: Modify C,^ by refining all elements such that rfg_^ ]\k > stol+ and 

coarsening all elements such that 7]g^ i\k < stol~. 

6 : if / ids > ttol'*^ then 

7: Ti ■<— ri/2. 

8 : end if 

9: Calculate 

10 : Calculate u\ from 

11 : end while 
12 : Calculate (5i. 

13: Multiply ttol’*', ttol“, stol’*', stol“ by the factor Gi. 

14: Setj= 0 , C' = C°. 

15: while ^j+i exists do 
16: j^j + 1. 

17: Tj+i =Ti. ^ 

18: Calculate from 

19: if / VT 2 ,j+i > ttol^ then 

Jp 

20 : Tj+i rj+i/ 2 . 

21 : Calculate from 

22 : end if 

|■p+'^ 

23: if /, ds < ttol then 

Jp 

24: Tj'+i ^ 2rj_|_i. 

25: Calculate from 

26: end if 

27: Create from C,^ by rehning all elements such that 'ffsxj+i\^ ^ stol'^ 

and coarsening all elements such that ? 7 |^< stol“. 

28: Calculate from 

29: Calculate 5j+i. 

30: Multiply ttol"*", ttol~, stol+, stol“ by the factor Gj+i. 

31: end while 

32: Output: j, t\ \\uh{P)\\^^(^^y 
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L°o(L°o) norm of the IMEX dG solution. 


5.5 Numerical Experiments 


We shall numerically investigate the presented a posteriori bound and the 
performance of the adaptive algorithm through an implementation based on the 
deal. II hnite element library m- All the numerical experiments have been 
performed using the high performance computing facility ALICE at the University 
of Leicester. For all the numerical experiments, we use polynomials of degree hve. 
Finally, we set ttol“ = 0.01 * ttol’*' and stol“ = 10“® * stol’*' as our temporal 
and spatial coarsening parameters. 


Remark 5.3. All unknown constants in the error estimators are set equal to one 
as is standard in a posteriori error analysis. It is believed that this is reasonable 


despite the fact that condition (5.21) is technically a strict limitation on whether 
or not we can continue our computations. 


5.5.1 Example 1 

Let G = (—4,4)^, e: = 1, a = (0,0)^, /o = 0 and uq = 106“^^^^+*^^^. The initial 
condition uq is chosen to be a Gaussian blob centred on the origin that is chosen 
‘large enough’ so that the solution exhibits blow-up; the blow-up set consists of 
a single point corresponding to the centre of the Gaussian. In order to observe 
how the error estimator behaves asymptotically, we choose a very small spatial 
threshold so that the spatial contribution to the error and the estimator are small. 
We then reduce the temporal threshold and see how far we can advance towards 


the blow-up time. The results are given in Table 5.1 


We know that asymptotically the solutions to (5.2) behave the same temporally 
as the solutions to (4.1) (at least in the case of zero convection) [6l]. This means 


that if our error estimator is good, we would expect to observe similar rates for A 
to those seen in Chapter 4. Although we do not know the blow-up time for this 


problem, we observe from Table 5.1 that 


||Wh||L°°(0,T;L°°(O)) CC 
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Table 5.1: Example 1 Results 


ttol+ 

Time Steps 

Estimator 

Final Time 

\\uh{T)\\L°°(n) 

1 

3 

9.5 

0.09375 

12.244 

0.125 

8 

24.6 

0.12500 

14.742 

(0.125)2 

19 

54.0 

0.14844 

18.556 

(0.125)3 

42 

66.7 

0.16406 

23.468 

(0.125)^ 

92 

218.5 

0.17969 

32.108 

(0.125)3 

195 

1142.4 

0.19043 

44.217 

(0.125)3 

405 

1506.0 

0.19775 

60.493 

(0.125)^ 

832 

1754.1 

0.20313 

83.315 

(0.125)8 

1698 

5554.2 

0.20728 

117.780 

(0.125)9 

3443 

6020.4 

0.21014 

165.833 

(0.125)3° 

6956 

33426.7 

0.21228 

238.705 

(0.125)33 

14008 

36375.0 

0.21375 

343.078 

(0.125)32 

28151 

66012.8 

0.21478 

496.885 

(0.125)33 

56489 

157300.0 

0.21549 

722.884 


From [M] , we know the relationship between the magnitude of the exact solution 
in the norm and the distance from the blow-up time. Thus, under the 

assumption that the numerical solution is scaling like the exact solution we get 

A(ttol''', iV) ~ ll'“/*llL°°( 0 ,r;L°°(O))- 

Therefore, we conjecture that 


A(ttol+,iV) oc 

This is obviously slower than the comparable results in Chapter 4 and a possible 
explanation for this will be discussed in the conclusions section. 

5.5.2 Example 2 

Let = (—4,4)^, e = 1, a = (1,1)^, /q = —1 and Uq = 0. This numerical 
example is interesting to study as not much is known about blow-up problems 
that contain convection. Here, the solution behaves as the solution to a standard 
convection-diffusion problem early on. As time progresses, the nonlinear term 
takes over and the solution begins to exhibit blow-up. As in Example 1, we 
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choose to use a small spatial threshold so the spatial contribution to the error and 


the estimator are negligible. We then reduce the temporal threshold and see how 


far we can advance towards the blow-up time. The results are given in Table 


5.2 


Table 5.2: Example 2 Results 


ttol+ 

Time Steps 

Estimator 

Final Time 

\\uh{T)\\L°°{n) 

1 

4 

3.6 

0.78125 

0.886 

0.125 

10 

3.6 

0.97656 

1.322 

(0.125)2 

54 

22.0 

1.31836 

3.269 

(0.125)3 

119 

47.5 

1.41602 

5.107 

(0.125)^ 

252 

132.1 

1.48163 

8.059 

(0.125)3 

520 

218.4 

1.51711 

11.819 

(0.125)3 

1064 

664.6 

1.54467 

18.139 

(0.125f 

2158 

1466.1 

1.56224 

27.405 

(0.125)8 

4354 

1421.7 

1.57402 

41.374 

(0.125)^ 

8792 

11423.0 

1.58243 

64.450 

(0.125)^° 

17713 

21497.8 

1.58770 

99.190 

(0.125)3^ 

35580 

21097.1 

1.59092 

145.785 

(0.125)12 

71352 

35862.0 

1.59299 

211.278 


From Table 5.2, we draw the conclusion that 


||Wh||L°°(0,T;L°°(O)) OC 


Although not much is known about blow-up problems that contain convection, it 
is reasonable to assume that because the nonlinear term dominates close to the 
blow-up time, the same relationship between the magnitude of the exact solution 
in the norm and distance from the blow-up time exists as in Example 1. 

If this is true, then under the same reasoning as in Example 1 we conclude that 


A(ttol+,iV) OC 


5.5.3 Example 3 

Let n = (—8,8)^, £ = 1, a = (0, 0)^, /o = 0 and the ‘volcano’ type initial condition 
be given by Uq = 10 -|- |/ 2 ^g-o. 5 (a: 2 +j/ 2 )^ blow-up set for this example is a 

circle centred on the origin — this induces layer type phenomena in the solution 
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around the blow-up set as the blow-up time is approached making this example 
a good test of the spatial capabilities of the adaptive algorithm. Once more, we 
choose a small spatial threshold so that the spatial contribution to the error and 
the estimator are negligible. We then reduce the temporal threshold and see how 


far we can advance towards the blow-up time. The results are given in Table 5.3 


Table 5.3: Example 3 Results 


ttoR^ 

Time Steps 

Estimator 

Final Time 

\\uh{T)\\L°°{n) 

8 

3 

15 

0.06250 

10.371 

1 

10 

63 

0.09375 

14.194 

0.125 

36 

211 

0.11979 

21.842 

(0.125)2 

86 

533 

0.13412 

31.446 

(0.125)3 

190 

971 

0.14388 

45.122 

(0.125)^ 

404 

1358 

0.15072 

64.907 

(0.125)3 

880 

5853 

0.15601 

98.048 

(0.125)3 

1853 

10654 

0.15942 

146.162 

(0.125)^ 

3831 

21301 

0.16176 

219.423 

(0.125)8 

7851 

143989 

0.16336 

332.849 

(0.125)9 

16137 

287420 

0.16442 

505.236 

(0.125)13 

32846 

331848 

0.16512 

769.652 

(0.125)11 

66442 

626522 

0.16558 

1175.21 


Once again, the data implies that 




Arguing as in Example 1, we again conclude that 


A(ttol+,iV) oc 


The numerical solution at t = 0 and t = T from the hnal numerical experiment 
(ttol+ = (0.125)^^) is shown in Figurethe corresponding meshes are displayed 
in Figure |5.2 The initial mesh has a relatively homogenous distribution of 
elements which is to be expected since the initial condition is relatively smooth. 
In the hnal mesh, elements have been added in the vicinity of the blow-up set and 
removed elsewhere, notably near the origin. The distribution of elements in the 
hnal mesh strongly indicates that the adaptive algorithm is adding and removing 
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Figure 5.2: Example 3: Initial (left) and final (right) meshes. 



Figure 5.3: Example 3: Initial (left) and hnal (right) solution prohles. 


elements in an efficient manner. 


5.6 Conclusions 


The error estimator produced performed adequately in both numerical examples 
but the blow-up time was approached at a much slower rate than expected given 
the results in Chapter 4. The reason for this lies in the signihcant differences 


between the 6 equations (4.12) and (5.22) which can in turn be traced back to the 


error equations (4.5) and (5.15). For a quadratic nonlinearity, the highest order 


error term in both of these equations is one power higher than the remainder 
of the error terms which can all be explicitly bounded by Gronwall’s inequality 
(heuristically, the ODE and PDE error analysis are still ‘equivalent’ at this point). 
However, in the PDE analysis we have no way of dealing outright with an norm 
necessitating the use of the Gagliardo-Nirenberg inequality. After application of 
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Gronwall’s inequality, a term of the form 


exp I J I iVed I ds 


(5.26) 


remains. With no way to estimate this directly through a continuation argument, 
we are forced to use the Cauchy-Schwarz inequality which destroys a large amount 
of information (we do something slightly different to get a better 6 equation but 
the end result is still the same — a loss of information caused by the necessity 
of having norms that are compatible with our continuation argument). If all the 
norms are the same (as in Chapter 4), there is no loss of information. Therefore, 
we conjecture that conducting an error analysis for the norm may lead 

to a recovery of the rates seen in Chapter 4. 
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Chapter 6 

A posteriori error estimation and 
adaptivity for a class of nonlinear 
interface problems 


6.1 Model problem 


We shall consider a model problem that is a simplihcation of the models given in 
|251 12B] which were constrncted to model the mass transfer of solntes throngh a 
semi-permeable membrane. The derivation of the models in [2S112S], based npon 
the works [SU [73l El 110011114] , shall not be restated here becanse we are primarily 
interested in the mathematical model. 

The compntational domain is snbdivided into two snbdomains fli and 
snch that SI = U T* where := f2i U f22 and Tj = dQ\dfl is the interface 
between the two snbdomains. To simplify things, we assnme that the interface is 
a non-intersecting piecewise linear cnrve. 

For T > 0, consider the model problem of hnding m : x (0, T] —)• M snch that 

du 

— - eAu + SL-'Vu + f(u)=0 inf2x(0,T], 

dt ^ ^ ( 6 . 1 ) 

u{-, 0) = uo in 


We need to angment (6.1) with snitable bonndary conditions. To that end, we 
split the bonndary dfl = TU Twhere T^i is the Dirichlet bonndary and T n 
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is the Neumann boundary. It is assumed that the intersection of the Dirichlet 
boundary with each of the subdivision boundaries has positive one-dimensional 
Hausdorff measure. We then impose boundary conditions for all t G (0,T]: 


u = 0 
eVu ■ Yi = g 
{eVu — aw) ■ n = 


on r^), 

on Tat n dflout, 
on Ttv n dQin- 


( 6 . 2 ) 


In addition to boundary conditions, we need to augment (6.1) with interface 
conditions. To that end, we require that the following equalities are satished 
across the interface for all t G (0,T]: 


(eVu - aw) ■ n|oj = p(m|q 2 - ulnj - + tC 2 M|o 2 )(a ■ n)|oi, 

(6.3) 

{eVu - au) ■ n|o 2 = p{u\n^ - u\n^) - r{wiu\Q^ + W 2 u\q^){sl ■ n)|n 2 , 

where p > 0 is the permeability coefficient, r G [0,1] is the friction coefficient and 
wi,W 2 are weights that satisfy wi + W 2 = I and 0 < < 1- The weak form 

of the model then reads as follows: find u G L^(0, T; fl (0,T; 
such that for almost every t G (0,T] we have 

+ B{u,v) + {f{u),v) =l{v) yv e (6.4) 

with 


B{u, v) := / (eVu — aw) ■ Vu dx — jV ■ auv dx + 

Jn Jn 

+ / p[u\-[v]ds+ / r{u}u,[av] ds, 


a ■ nuvds 




(6.6) 


l{v) := / gvds, 
Jvm 


where {u}w ■= wiu\n-^ -|- W 2 u\n 2 is the weighted average of u across Tj. We make 
the following assumptions upon the data: g G L^(rAr), Uq G 0 < e < 1, 
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a G ^ and / : M —)■ M must satisfy the growth condition 


\f{u) — f{v)\ < L\u — v\{l + \u\ + \v\)'^ Wu, n e M, 


( 6 . 6 ) 


for some L > 0 and 0 < fi < 2. We hnish the section by introducing a coercivity 
result for the bilinear form B. 


Theorem 6.1. Let c* denote the constant in the trace inequality (Theorem 2.10) 
and define Ai := ||a||£cx>(r;). For any v G the bilinear form B satisfies 

B{v,v) > ^|||n||p + essjnf(-V ■ a) - c^arwAi{l + 4:C^arwAie~^)^ ll^lP, 


where a™ := 2 — W 2 I + max 


rwi - - 


rw2 - - 


and 


l^^lll := ( ^l|Vn|| + - 


|a-n|n^ds + J p\[v]\‘^ds\ 


Proof. The multidimensional integration by parts formula implies that 

I {siv) ■ Vv dx + - i V ■ slv^ dx = - j a.-nv'^ds+ j {v}[siv]ds. 

Jn ^ Jn 2 J-p. 

Application of this to B{v,v) yields 

> |||n|||^ + ^essinf(—V ■ a)||r>||^ + j {r{v}w — {v})[av]ds. 

2 o Jp. 

Denoting v evaluated on Di by Vi and v evaluated on D 2 by V 2 then using Young’s 
inequality we obtain 


- {i^})[an] ds 


< 


rwi 


j |a ■ n|r;^ ds + rw 2 — - j 

OTi 2 Jp 


|a ■ n|n 2 ds 


+ r\wi-W 2 \ / |a ■ n||r;i||n 2 | ds 


< arwAi ( I vf ds + J V 2 ds 
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Using the trace inequality with S = -c^ we obtain 


- {^^})[an] ds 


< c*ar^^i(5||Vn||^ + (l + 5 ^)||'y||^) 

^ tIII'^IIP + + 4 :C^arwAi£~^] 


IniP. 


Therefore, combining the results 

B{v,v) > |||n||p + ^essjnf(—V ■ a)||r>|p — 




{r{v}^ - {n}) [an] ds 


> -lllnl 


+ ess^inf(—V ■ a) — c^arwAi[l + 


This completes the proof. 




□ 


6.2 Space-time discretisation 

The discontinuous nature of the model on the interface makes a spatial dG 
discretisation a natural choice. Thus, with the same reasoning as in the previous 


chapter, we consider an IMEX dG discretisation of problem (6.4). 


To that end, consider a subdivision of [0,T] into time intervals of lengths 


U, 


Tr, 


such that = T for some n > 1 and set := 0 and := ■ 


i=i 


i=i 


Denote an initial triangulation by and associate to each time step A; > 0 a 
triangulation which is assumed to have been obtained from by locally 
rehning and coarsening All of the meshes are assumed to be aligned with 

the interface in the sense that no part of the interface is contained in the interior 
of any element and aligned with the boundary in the sense that the points of 
intersection between the Dirichlet and Neumann boundaries, if they exist, must 
all be at the vertex of an element. To each mesh we assign the finite element 
space := Vh{C^) given by (2.5) and we set := /(m^) for brevity. For each 
mesh let Vint denote the union of all interior edges that do not lie on the 

interface. Finally, for t G we let T denote the union of all edges in the 

mesh U that do not lie on the interface or Neumann boundary. 
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The IMEX dG method then reads as follows. Set to be a projection of uq 
onto 14 °. For /c = 1, n, hnd G Vjf such that 


fe-i 


Tk 


for all G Vjf where 


,vn +B{ulvt)+Mul,vt) + (6.7) 


■= -Vv^dx- ^ f V ■ Au^v^dx 

Ke(;kdK 


X 


EcrourM 


llL.k 

I E^E 


■ bS ds+ 


u, 


K&C 


k JdKout\Ti 


SLV^ ds 


+ [ pH] ■ H]ds+ [ r{ul}^[a.v'(\ds, 

Jr, Jr, 

K,H,;Q :=- Y. f ■ K] + ■ [“!:] <'*>■ 

T 7 ^T^\ \T. . ^ E 


ECTD^Tir^t ' 


(6.8) 


We shall take to be the orthogonal projection of Uq onto although other 
projections onto Vj^ can also be used. 

For /c > 0, the residual Rk is dehned on the interior and edge of an element 
77 G as follows 


—^^- 4 *-^ ~ ^ 

in K 

0 

on dK n Td 

g — eVu^ ■ n 

on (977 n (Tat n dflout) 

g — {eVu^ — a.u’l) ■ n 

on (977 n (Tat n dflin) 

(aw^ - eVul) ■ UK + 

on (977 n Ti if 77 C iR 

-ra-n;^{4}^ 


(a^ - eVul) ■ UK + p«|oi - u’^hlH 

on (977 n r* if 77 C ^2 

-ra-nif{4}^ 


-|[v4] 

on (977 n Tint 
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6.3 An a posteriori bound for the IMEX dG method 

At each time step k, we decompose the dG solution into a conforming part 
c G and a non-conforming part G such that = u’l c + u^d- 

Given t G we dehne Uhit) to be the linear interpolant with respect to t 

of the values and viz., 

Uhit) := lk-iit)ul~^ + lkii)ul. 


We dehne Uh,cit) and Uh,dit) analogously. We can then decompose the error 
e := u - Uh = ec - Uh,d where Cc ■-= u - Uh,c- 

Lemma 6.1. Given t G then for any v G Hjjikl) we have 


+Bie,v) + ifiu) - fiuh),v) = liv) - + fiuh),v ) - Biuh,v). 


Proof. Follows from (6.4). 


□ 


From Lemma 6.1, it follows that 


de \ 

— , v \ + S(e, v) + (fiu) - fiuh),v) = - fiuh),v) + B{ul, v) 

- B{uh,v) +l{v) - + -B{ul,v). 


(6.9) 


Finally, we use (6.7) to conclude that for any G 

(1^ w) + B{e, v) + {f{u) - fiuh),v) = - fiuh),v) + B{uh, v) 

- B{uh,v) +l{v-v’^) - + f-\v - B{ul,v - vl) 

+ Kh{ulA)- 


We are now ready to state our a posteriori estimator. Due to the nature of the 
error bound to be presented, it is easier to separate the estimator into two parts. 
As in previous chapters, a subscript S denotes parts of the estimator related to 
estimating space while a subscript T denotes parts of the estimator related to 
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estimating time. In this way, for t G (t^ , i ]a is given by 


where 


with a 
for t e 


VA ■= VSi,k + VS2,k + VS3,k + VS4,k + hTi,fc + VT2,k + VT3,k, 


VSi,k ■ = ( ^\\^k\\l 2^K) + ^ W^kWh^E) 

K&C,^ E<ZdK 




E EIIK]'" 


hp, w r fcl I |2 


1/2 


EcroUFi^t 


h 


E 'I-E r ^ 

— [aw; 


\L^E) ^ ^ 

ECFm 

hr 


h\ I \l‘^{E) 
\ 1/2 


vs2,k ■■= (5: ^IIK]|li.(i,) + E -^\\[^u,]\\l2^E) ) , 


kEcT 


EdV 


J 

1/2 


^ ‘^Af|IMIlL2(E) 

^EcTNf^dnout KdE Ecitfinr 

1/2 

EE E “piiMiit(E)' 

^EcFi XcEEcifenr 

VTi,k ■= I - Uh) - - “/i) 11 > 

2 

|a ■ n| — Uh\ ds \ , 


hSs.fc • = 


h54,fc : = 


VT2,k ■ = 


'rMndiiout 




dT3,k ■■= +^P ^^‘^\^\\{ul-Uh] 


w\ WL^iTi) ’ 


, := 2p + 2r^p ^ m&yi{wl,wl]A^i and An ■= ||a||i,oo(rj^n 90 o„t)- Similarly, 
pb is given by 


Vb ■— VS5,k + VSe,k + VT4,k, 
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where 


1/2 


VS5,k-= 1 ^ ^ ^KhE\\[Uh]\\l2i^E) 

^xec'^-^uc'' ec-R-b 

..fe-n 2 \V2 


^Ecr 




T-fc 


E2(E), 


hS6,fc - 

m.t ■■= 11/‘“‘ - f{uh) - V ■ a(u‘ - m) 11, 


with 


(Ti^ := max {L,2'^ (l + 2| |m/i| |ioo(^) + 


L°°(EBnr) 


)^ + 11 V ■ a| 


The hrst three terms on the right-hand side of (6.10) approximate the temporal 


part of the error. To begin bounding, we decompose these terms, viz., 

(/"■' - f{uf,),v) + B{ul v) - B{uh, v)=T^+T^ + T^ + T 4 , ( 6 . 11 ) 


where 


Ti := ^ / {eS/{ul-Uh) -&{ul^-Uh))-Vvdx, 




T2-.= Y1 / (/'■'-/K)-V-a(4-w,))ndx, 


iC'GC^-^UC^ 


T 3 := / a ■ n(M^ - M/i)r>cis, 

-'Tjvnen oiif 

^4 := / - Mh] ■ [n] ds[an] ds. 

dr, dPi 

Bounding Ti requires a simple application of the Cauchy-Schwarz inequality: 


lEI < 


( 6 . 12 ) 


T 2 is also bounded by the Cauchy-Schwarz inequality, viz.. 










Ts is bounded by the Cauchy-Schwarz inequality as follows 


iTsl < riT2,k ( [ |a-n|u^ds^ 

\J r j^ndfiout / 

^ vt^AMW- 


(6.14) 


Finally, T 4 is bounded by the Cauchy-Schwarz inequality, viz.. 


ITaI < vtsAWA 


(6.15) 


The remainder of the terms on the right-hand side of (6.10) give rise to parts 


of the space estimator. We start by noting that, through application of the 
multidimensional integration by parts formula on each element, we have 


{v - vl) - f \v- - B(uA u - uQ = T 5 + Tg, (6.16) 


where 


T, : = 


f Rk{v - v’A) dx + Y^ f Rk{v-v’A) dx, 

EcdK 


^6 •= Y1 / - ^h) ds. 


Ecr^„t' 


We set G fl to be the hnite element interpolant from Theorem 


2.2 Application of the Cauchy-Schwarz inequality together with the interpolation 
estimates from Theorem |2.2| yields 


hlr 


1/2 


1/2 


ir^l < 5: -fWR 






hr 




k\\L^(E) 


E^dK 


A, VSi,k 


E 

1/2 


1/2 


k E 


1/2 


£ 

He 


, KeC*’ EcdK 


u \\2 

L^lE) 


(6.17) 


1/2' 


,icec'= 


2 

L^K) 


Y Y 

K€(>‘ EcdK 


2 

LHE) 


< 


hSi,fc||T| 
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Tg is also bounded through the interpolation estimates of Theorem 2^ together 
with the Cauchy-Schwarz inequality, viz., 


\n\ < 


< 

rs_/ 


< 


KecTm / \Ecri, 

VSi,k( 

\Ecr^^t J 


1/2 


^11 1-1 |2 


L^E) 


(6.18) 


vsiMW^l 


Finally, Kh{u^,v^) is bounded through the Cauchy-Schwarz inequality and the 
inverse estimate along with the shape-regularity of the mesh as follows 

|/f4<,«ai< 5^ f £|V«‘||[4]|£is 

EcrourM 


< 


< 


< 


\£;crDuri„t ® / \Ecrouri„t / 

VEcrouri^t / 


(6.19) 


VSi,k 


hSi,fc||Tl 


Putting together all these results we obtain 


f)p \ 

+B{ec,v) < \{f{u) - f{uh),v)\ + 


du 


h,d 


dt 




+ \B{uh,d,v)\ 


( 6 . 20 ) 


+ {VSi,k + VTr,k + VT 2 ,k + VT3,k)\\\v\\ \ +Vn,kM 


To bound B{uh,d,v), we note that 


B{uh,d, v) — T'l + Tg + Tg + Tio, 


(6.21) 






where 


Tv := 


T« := 


/ {eVuh,d - ^Uh,d) ■ Vv dx, 
/ W ■ aiUh,dvdx, 

rk-ii ,rk 


-R'sC'^-^uC' 


To := 




a ■ YiUh^v ds, 


^gj^fe-iy^fe •^s^^outnrjv 

Tio := / p[Mh,d] ■ H (is+ / hs. 

Jri Jvi 


To bound T 7 , we use the Cauchy-Schwarz inequality along with the estimates from 
Theorem 12.31 to conclude that 


I^tI < vs2,k\\M\, 


( 6 . 22 ) 


while Tg is bounded through Holder’s inequality, the Cauchy-Schwarz inequality 
viz.. 


and Theorem 2.3 


1/2 

iTgl < I 11^ ■ 

^ hss.fcibll- 


(6.23) 


Tg is bounded using Holder’s inequality, the Cauchy-Schwarz inequality, the trace 
inequality (with 6 = hx) and the bounds from Theorem |2.3| along the shape-regularity 
of the mesh as follows 


\T9\< 


/ \a ■ n\\uh,d\‘^ ds 

J dKoutn^ff 

\ 1/2 

^ Il2(e) I iii'^iii 

EcrMndiiout / 

I y~^ I 17^2(7^') -|- \ Wh,d\\L'^{K') 

.EcTMf^dUout KcE 


(6.24) 


1/2 


< 


hS3,fc|lTl 






Finally, Tio is bounded using Holder’s inequality, Young’s inequality, the 
Cauchy-Schwarz inequality, the trace inequality and the bounds from Theorem 


2.3 along the shape-regularity of the mesh, viz.. 


iT’iol < p\[uh,d\\‘^ ds + j r‘^p ^A^i\{uh,d}w?d^ |||n| 

< + ^ 1 ^ 2 ) ds^ 11 |n| 11 


yEcFi kcE 


1/2 


(6.25) 


< 


hS4,fc|lTl 


To bound the remaining nonconforming term, we use the Cauchy-Schwarz inequality 
and the bounds from Theorem 12.31 as follows 


du 


h^d 


dt 


< 


du 


h,d 


dt 


^^11 ^ VSe,kM 


(6.26) 


Combining these results, using the dehnition of our estimators and setting v = Cc 
we obtain 

^^(lledP) +B{ec,ec) < |(/(m) - f{uh),ec)\+PA\\\ec\\\ +PB\\ec\\- (6-27) 


We must now bound the nonlinear term. The growth condition (6.6) and the 
triangle inequality imply that 


|(/(m) - f{uh),ec)\ <L |e||ec|(l -f |m| \uh\)'" dx 


(6.28) 


<L (led |Mh,d|)|ec|(l + 2\uh\ + \uh,d\ + \ec\Y dx. 
Jn 


Thus, using the power mean inequality we have 


I (/('^) ~ 5 ^c) I Y Til + Ti2 + Ti3 -|- Ti4, 


(6.29) 
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where 



Ti 2 :=max{L, 2 '" ^L] f \uh,d\\ec\^^'^ dx, 

Jn 

Ti 3 := max {L, 2^~'^L} /(I + 2\uh\ + \uh,d\Y\ec\'^ dx, 

Jo, 

Ti4 := 

To bound Tn, we use the Cauchy-Schwarz inequality, Holder’s inequality and 
Theorem 12.31 to conclude that 


^11 ~ hSs.fclkcll, 


(6.30) 


while Ti 2 is bounded through Holder’s inequality, Theorem 2^ and U’ embeddings 


VVAJ.XXV.. -*-1^ L..»V-/CXXXVJ.V_.VJ. UXXXV^l.X^XX J.J.V^XVJ.V.X O XXXV..V^l.XCXXXOJ , J. J 

if 0 < /r < 1 or the Gagliardo-Nirenberg inequality 


if 1 < /i < 2 , viz.. 


Ti 2 <max{L, 2 ^ ^L}\\[uh]\\L°°{T)\\ec\\^^^ if 0 </i < 1 , 


ri 2 <max{L, 2 ^ ^L}\\[uh]\\L--{r)\\ec\f\\Vec\\'" ^ if 1 </i < 2 . 
To bound T 13 , we use Holder’s inequality and Theorem 2^ 


as follows 


/ill I |2 

L«)(r)l NCr 


Ti3<max|L,2^ ^L} (l + 2| |m/i| + | |l°°(o)) II^cII^ 

< max{L, 2 ^"^L}(l + 2\\uh\\L^(n) + IIMIlL-(r)) 

Finally, the Gagliardo-Nirenberg inequality implies that 

Ti4<max{L,2^-iL}||e,|n|Vee| 


(6.31) 


(6.32) 


(6.33) 


Let C and K denote generic positive constants and dehne ai '■= max{2L,2^L} 


for brevity. Applying the above bounds to (6.27) and using the coercivity of the 
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6.1) and Young’s inequality yields 


bilinear form B (Theorem 


4(lkcin + IlledlP < + Cr^slledl + crilledl^’^''+ (dn + crallVedl^ ^ 

ut (6.34) 


where 


+ KaL\\Ve,\n\\e,\\^ 


CTi := 


(T2 := 


KaL\\[uh\\\L°°{r) if 0 </i < 1 

0 if 1 < /i < 2 ’ 

0 if 0 < /i < 1 

KaL\\[uh\\\L°-(r) ifl</i<2’ 

an := aL(l + 2\\uh\\L°°iQ) + -^11MIli”(r))^ - essjnf(-V ■ a) 

+ 2c*Q;ru)*^i (l + 

Another application of Young’s inequality yields 

) + |||ec|||^ < C[r]\ + T?]!) + ailledd'^'^ + + an 

+ CX 2 II Veel r-'+ ifail I Veel I led d- 


dr"®"' 


(6.35) 


Application of Gronwall’s inequality to (6.35) together with the bound 


l|ec(0)|d< ||e(0)|d+esssup||Mft,dr r l|e(0)|d+esssup V hi?||[Mdlli 2 (ij), (6.36) 

0<t<T 0<t<T 

implies that for any t G [0, T] we have 

I |ec(t) 11^ < CH{t)G (^0 + I led ds^ , (6.37) 


where || dl* is the + L°°{L‘^) type norm 

ir(^)ll* •= ( Inli“(o,t;L2(n)) + 


pt \ 1/2 

I |||M||d(isJ , 
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and 


H{t):=exp^J cr2||Vec||^ ^ ds + Kai j ||Vec||^ds^, 

G := exp ^ y (Jndsy 

/•T /.T 

(/) := ||e(0)|p + / r]\ds + T 77 I ds + esssup ^/ie||[m/j]||^ 2 (e)- 
7o Jo 0 <t<T 


In order to construct a practical error estimator from (6.37), we employ a continuation 
argument. To that end, we define the set 


I:={te[0,T] 1 \\eM\l<SG<p}, 


where d > G is a parameter to be chosen. Clearly I is bounded; furthermore, we 
know I is non-empty because 0 G /. Let t* denote the maximal value of t in / 
and assume that t* < T. We proceed as in previous chapters by bounding the 


remaining error terms in (6.37). Firstly, Holder’s inequality implies that 

[ < ||ec(t*)||i’''^ f cTids<(dG0)^ f aids, (6.38) 

Jo Jo Jo 

while through the Cauchy-Schwarz inequality and embeddings we obtain that 

rT \ 1/2 / X 1/2 


(^2 


iVedli* ^ ds < a^ds^ UVedP^ ^ ds^ 


< T^-2 


< T^-2 


rT \ 1/2 


M-i 

2 


(6.39) 


a, 


ds] (dG0)V. 


'0 / 

Finally, we use the properties of embeddings to conclude that 

\ G2 


I Verl li* ds < T 


/ rr \ 

^~Hj llVedpdsJ <T^-^5G(|))G\ 


(6.40) 
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Putting these results into (6.37), we conclude that 


1 +M 


where 


\e,it*)\\t<C^/JG[(j)+{dG(j))^ aids), 


ij :=exp alds^ ' (dG(/.)V + (5G0)/^/2 j . (6.42) 


(6.41) 


Now, suppose that the upper bound in (6.41) is strictly less than the upper bound 
of the set /, viz.. 


G'lpci^ + {dG(j)) G j aids ] < SGcf), 


1+M 


or equivalently. 


1+M 


Gijj ( (j) + {dG(j)) 2 ai ds ] < 6(j), 


(6.43) 


(6.44) 


then t* cannot be the maximal value of t in / because we just showed that | |ec(t*) |\l 
satishes a bound strictly less than that assumed in the set / — a contradiction. 


Therefore, providing (6.44) is satished, I = [0,T] and we have our desired error 
bound once we select 6. Taking the limit, we can select S to be the minimiser of 


1 +^ 


G^lJ[(j)+{dG(j))^ j aids] -6(j) = 0, 6>G. 


(6.45) 


In order to state the hnal theorem, we need to extend the energy norm to include 
functions in I 4 . To that end, for t G we (re)de£ne 




( ^\\^^\\l\K) + \ [ |a-n| 

\A'eC'=-iuc'“ 


n^ds + J p|[n]|^ds 




\ 1/2 


Ecr ^ Ecr 

We then have the following result. 
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Theorem 6.2. The error of the IMEX dG discretisation of problem (6.4) satisfies 


k(r)||. < 


provided that the solution to (6.45) exists. 


Proof. From the triangle inequality, we have 

||e(T)||, < ||e,(r)||, + \MT)\U < + \MT)\U. 

Thus, all that remains is to bound \ \uh,d(T)\\^,; the part of this term was 

bounded in (6.36) while the part of the term was bounded in (6.22), (6.24) 


and (6.25). Thus, 


IKd(r)ll. < 


This completes the proof. 


□ 


6.4 Numerical experiments 

We shall numerically investigate the presented a posteriori bound through an 
implementation based on the deal. II hnite element library m- In particular, 
we shall use Algorithm 3.1 from Chapter 3. Spatial refinement and coarsening 
are driven by the term subject to a spatial rehnement threshold stol+ and 
a spatial coarsening threshold stol“. As in Chapter 3, we define 


VT,k ■= / {vn,k + r]T2,k + Vn,kYdt + T ril dt, 

Jtk-i Jtk-i 


(6.46) 


the sum of which bounds the full time estimator. Temporal rehnement is then 
carried out using 7)T,k subject to a temporal threshold ttol on each time interval. 

In all our numerical experiments, we use polynomials of degree two and an 
initial 4x4 uniform quadrilateral mesh. Finally, the spatial coarsening threshold 
is set to stol“ = 0.001 * stol'*’. 
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Figure 6.1: Example 1: Meshes produced by the adaptive algorithm for e = 0.1 
(left) and £ = 10“^ (right). 



Figure 6.2: Example 1: Solution prohles for e: = 0.1 (left) and £ = 10 ^ (right). 



6.4.1 Example 1 

Let = (-1,0) X (-1,1), = (0,1) X (-1,1), a = (1,1)^, f = -1, uo = 0 

and T = 1. For the interface parameters, we set p = 0.1, r = 0.5, wi = 1 and 
tC 2 = 0. Under this choice of interface parameters, the solution to (6.4) exhibits 
both boundary and interface layers of width 0{e). Solution prohles and meshes 
produced by the adaptive algorithm at the hnal time are given in Figures 6.1 and 


6.2, respectively. The meshes generated by the adaptive algorithm clearly show 


that the error estimator is picking up both the interface layer and the boundary 
layer. 

To observe the rates of convergence of the error estimator 0, we begin by hxing 
a small temporal threshold; the spatial threshold is then reduced to observe the 
spatial rates of the estimator. We then £x a small spatial threshold so that all 
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Figure 6.3: Example 1: Spatial and temporal rates. 


layers are sufficiently resolved and reduce the temporal threshold to observe the 


temporal rates of the estimator. The results, given in Figure |6.3[ show that the 
space and time estimators are of optimal order. 


6.5 Conclusions 


We derived an a posteriori error estimator for a nonlinear interface problem that 
is used to model the flow of solutes through semi-permiable membranes. The 
error estimator displayed optimal spatial and temporal rates under Algorithm 
3.1. Furthermore, the error estimator was able to detect and rehne the interface 
layer without wasting degrees of freedom on the opposite side of the interface. The 
constant from Gronwall’s inequality is of the order exp(£“^), which is impractical 
for the convection-dominated regime. A different treatment of the interface terms 


in Theorem |6.1| may yield a tighter error bound, but it is not currently clear how 
to address this issue. 
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Chapter 7 


Summary and outlook 


The aim of this work was to advance the understanding of adaptive algorithms 
for spatial hnite element discretisations of parabolic problems — this was achieved 
in two ways. Firstly, in Chapter 3 an adaptive algorithm was proposed that 
utilised an error estimator derived for a backward Euler dG discretisation of a 
linear non-stationary convection-diffusion equation. This adaptive algorithm was 
applied to test problems in Chapter 3 as well as to a nonlinear interface problem in 
Chapter 6; in all test cases the error estimators were reduced at the theoretically 
expected rate with respect to the discretisation parameters. Secondly, adaptive 
algorithms designed to converge to the blow-up time of an ODE with polynomial 
nonlinearity were explored in Chapter 4. This led to the development of an 
adaptive algorithm in Chapter 5 that was designed to approximate the blow-up 
time of a semilinear parabolic PDE. The adaptive algorithm was then applied in 
two numerical experiments and shown to approximate the blow-up time of both 
problems. We shall now discuss some ways in which the results of this work could 
be extended on a chapter by chapter basis. 

In Chapter 3, we derived an a posteriori error estimator for a backward Euler 
dC discretisation of a linear non-stationary convection-diffusion equation and 
developed an adaptive algorithm to utilise the error estimator. There are several 
ways that the work in this chapter could be extended: 

• The extension of the error estimator to include a variable diffusion coefficient. 

• The extension of the error estimator to higher order time stepping schemes. 



• The proof of lower bounds for the given a posteriori error estimator. 

• A rigorous proof that the adaptive algorithm minimises the spatial and 
temporal parts of the estimator. 

In Chapter 4 and Chapter 5, we investigated the numerical approximation 
of blow-up through a posteriori error estimation and looked into the creation of 
adaptive algorithms designed to approximate the blow-up time. The work in these 
chapters could be furthered by: 

• The extension of the error estimators to include more general nonlinearities. 

• The extension of the error estimators to higher order time stepping schemes. 
In particular, it would be of great interest to study hp time stepping schemes 
for blow-up problems. 

• Conducting the error analysis for a different norm. In particular, conducting 
an error analysis for the L°°{L°°) norm may yield a faster approach to the 
blow-up time. 

Finally, in Chapter 6 we derived an a posteriori error estimator for an IMEX 
dC discretisation of a nonlinear interface problem. The error estimator was then 
applied to a test problem using the adaptive algorithm from Chapter 3. The work 
in this chapter could be extended by: 

• Removing or weakening the exponential dependence on £ from the error 
estimator, possibly via a spectral estimate. 

• The extension of the error estimator to include variable diffusion, time 
dependent coefficients and data that is (possibly) discontinuous across the 
interface. 

• The extension of the error estimator to the full system of equations considered 

in [251 [2S]- 


99 


Bibliography 


[1] Robert A. Adams and John J.F. Fonrnier. Sobolev spaces, volnme 140. 
Academic press, 2003. 

[2] Martial Agueh. Gagliardo-Nirenberg ineqnalities involving the gradient 
L^-norm. Comptes Rendus Mathematique, 346(13):757-762, 2008. 

[3] Mark Ainsworth and J. Tinsley Oden. A posteriori error estimation 
in finite element analysis. Pnre and Applied Mathematics (New York). 
Wiley-Interscience [John Wiley & Sons], New York, 2000. 

[4] Georgias D. Akrivis, V.A. Dongalis, Ohannes A. Karakashian, and W.R. 
McKinney. Nnmerical approximation of blow-np of radially symmetric 
solntions of the nonlinear Schrodinger equation. SIAM Journal on Scientific 
Computing, 25(1):186-212, 2003. 

[5] Rodolfo Araya, Edwin Behrens, and Rodolfo Rodriguez. An adaptive 
stabilized hnite element scheme for the advection-reaction-diffusion 
equation. Applied Numerical Mathematics, 54(3):491-503, 2005. 

[6] Rodolfo Araya, Abner H. Poza, and Ernst P. Stephan. A hierarchical 
a posteriori error estimate for an advection-diffusion-reaction problem. 
Mathematical Models and Methods in Applied Sciences, 15(07):1119-1139, 
2005. 

[7] Douglas N. Arnold. An interior penalty hnite element method 
with discontinuous elements. SIAM journal on numerical analysis, 
19(4):742-760, 1982. 


100 



[8] Louis A. Assale, Theodore K. Boni, and Diabate Nabongo. Numerical 
blow-up time for a semilinear parabolic equation with nonlinear boundary 
conditions. Journal of Applied Mathematics, 2008, 2009. 

[9] Garth A. Baker. Finite element methods for elliptic equations using 
nonconforming elements. Mathematics of Computation, 31(137):45-59, 
1977. 

[10] J.M. Ball. Remarks on blow-up and nonexistence theorems for nonlinear 
evolution equations. The Quarterly Journal of Mathematics, 28(4):473-486, 
1977. 

[11] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II—a general-purpose 
object-oriented hnite element library. ACM Trans. Math. Software, 
33(4):Art. 24, 27, 2007. 

[12] E. Bansch, F. Karakatsani, and Gh. Makridakis. The effect of mesh 
modihcation in time on the error control of fully discrete approximations 
for parabolic equations. Applied Numerical Mathematics, 67:35-63, 2013. 

[13] Soren Bartels. A posteriori error analysis for time-dependent 
Ginzburg-Landau type equations. Numerische Mathematik, 99(4):557-583, 
2005. 

[14] Soren Bartels and Riidiger Miiller. Quasi-optimal and robust a posteriori 

error estimates in for the approximation of Allen-Gahn equations 

past singularities. Mathematics of Computation, 80(274):761-780, 2011. 

[15] Amal Bergam, Ghristine Bernardi, and Zoubida Mghazli. A posteriori 
analysis of the hnite element discretization of some parabolic equations. 
Mathematics of computation, 74(251):1117~1138, 2005. 

[16] Marsha Berger and Robert V. Kohn. A rescaling algorithm for the numerical 
calculation of blowing-up solutions. Communications on pure and applied 
mathematics, 41(6):841-863, 1988. 


101 



[17] Stefano Berrone and Claudio Canute. Multilevel a posteriori error analysis 
for reaction-convection-diffusion problems. Applied numerical mathematics, 
50(3):371-394, 2004. 

[18] A. Bonito and R. Nochetto. Quasi-optimal convergence rate of an adaptive 
discontinuous Galerkin method. SIAM Journal on Numerical Analysis, 
48(2):734-771, 2010. 

[19] Alexander N. Brooks and Thomas J.R. Hughes. Streamline 
upwind/Petrov-Galerkin formulations for convection dominated flows 
with particular emphasis on the incompressible Navier-Stokes equations. 
Computer methods in applied mechanics and engineering, 32(1): 199-259, 
1982. 

[20] Chris J. Budd, Weizhang Huang, and Robert D. Russell. Moving mesh 
methods for problems with blow-up. SIAM Journal on Scientific Computing, 
17(2):305-327, 1996. 

[21] Z. Cai, X. Ye, and S. Zhang. Discontinuous Galerkin hnite element methods 
for interface problems: A priori and a posteriori error estimations. SIAM 
Journal on Numerical Analysis, 49(5):1761~1787, 2011. 

[22] Z. Cai and S. Zhang. Recovery-based error estimator for interface problems: 
Conforming linear elements. SIAM Journal on Numerical Analysis, 
47(3):2132-2156, 2009. 

[23] Zhiqiang Cai and Shun Zhang. Robust residual-and recovery-based a 
posteriori error estimators for interface problems with flux jumps. Numerical 
Methods for Partial Differential Eguations, 28(2):476-491, 2012. 

[24] A. Cangiani, E.H. Georgoulis, I. Kyza, and S. Metcalfe. Adaptivity and 
blow-up detection for nonlinear evolution problems. Preprint, 2015. 

[25] Andrea Cangiani, Emmanuil H. Georgoulis, and Max Jensen. Discontinuous 
Galerkin methods for mass transfer through semipermeable membranes. 
SIAM Journal on Numerical Analysis, 51(5):2911-2934, 2013. 


102 



[26] Andrea Cangiani, Emmanuil H. Georgoulis, and Max Jensen. Discontinuous 
Galerkin methods for fast reactive mass transfer through semi-permeable 
membranes. Applied Numerical Mathematics, 2014. 

[27] Andrea Gangiani, Emmanuil H. Georgoulis, and Stephen Metcalfe. Adaptive 
discontinuous Galerkin methods for nonstationary convection-diffusion 
problems. IMA Journal of Numerical Analysis, 2013. 

[28] J. Gascon, G. Kreuzer, R. Nochetto, and K. Siebert. Quasi-optimal 
convergence rate for an adaptive hnite element method. SIAM Journal 
on Numerical Analysis, 46(5):2524-2550, 2008. 

[29] J.H. Ghaudry, D. Estep, V. Ginting, J.N. Shadid, and S. Tavener. A 
posteriori error analysis of IMEX multi-step time integration methods for 
advection-diffusion-reaction equations. Submitted for publication, 2014. 

[30] Zhiming Ghen and Jia Feng. An adaptive hnite element algorithm with 
reliable and efficient error control for linear parabolic problems. Math. 
Comp., 73(247): 1167-1193 (electronic), 2004. 

[31] James Goleman and Gatherine Sulem. Numerical simulation of blow-up 
solutions of the vector nonlinear Schrodinger equation. Physical Review E, 
66(3):036701, 2002. 

[32] Javier De Frutos, Bosco Garcia-Archilla, and Julia Novo. A posteriori error 
estimates for fully discrete nonlinear parabolic problems. Computer methods 
in applied mechanics and engineering, 196(35):3462-3474, 2007. 

[33] Arturo de Pablo, Mayte Perez-Llanos, and Raiil Ferreira. Numerical blow-up 
for the p-Laplacian equation with a nonlinear source. In Proceedings of the 
11th International Conference on Differential Eguations (EguadiffOS), pages 
363-367, 2005. 

[34] A. Demlow and E. Georgoulis. Pointwise a posteriori error control for 
discontinuous Galerkin methods for elliptic problems. SIAM Journal on 
Numerical Analysis, 50(5):2159-2181, 2012. 


103 



[35] Stefka Dimova, Michael Kaschiev, Milena Koleva, and Daniela Vasileva. 
Numerical analysis of radially nonsymmetric blow-up solutions of a 
nonlinear parabolic problem. Journal of computational and applied 
mathematics, 97(l):81-97, 1998. 

[36] Vit Dolejsi, Alexandre Ern, and Martin Vohralik. A framework for robust a 
posteriori error control in unsteady nonlinear advection-diffusion problems. 
SIAM Journal on Numerical Analysis, 51(2):773-793, 2013. 

[37] W. Dorfler. A convergent adaptive algorithm for Poissons equation. SIAM 
Journal on Numerical Analysis, 33(3):1106-1124, 1996. 

[38] Sever Silvestru Dragomir. Some Gronwall type inequalities and applications. 
Nova Science Publishers, 2003. 

[39] Todd Dupont. Mesh modihcation for evolution equations. Math. Comp., 
39(159):85-107, 1982. 

[40] Kenneth Eriksson and Claes Johnson. Adaptive hnite element methods for 
parabolic problems I; A linear model problem. SIAM Journal on Numerical 
Analysis, 28(l):43-77, 1991. 

[41] Kenneth Eriksson and Claes Johnson. Adaptive hnite element methods for 

parabolic problems II: Optimal error estimates in and SIAM 

Journal on Numerical Analysis, 32(3):706-740, 1995. 

[42] Kenneth Eriksson and Claes Johnson. Adaptive hnite element methods for 
parabolic problems IV: Nonlinear problems. SIAM Journal on Numerical 
Analysis, 32(6):1729-1749, 1995. 

[43] Kenneth Eriksson and Claes Johnson. Adaptive hnite element methods for 
parabolic problems V: Long-time integration. SIAM journal on numerical 
analysis, 32(6):1750-1763, 1995. 

[44] Kenneth Eriksson, Claes Johnson, and Stig Larsson. Adaptive hnite element 
methods for parabolic problems VI: Analytic semigroups. SIAM journal on 
numerical analysis, 35(4):1315-1325, 1998. 


104 



[45] Alexandre Ern and Jennifer Proft. A posteriori discontinuous Galerkin error 
estimates for transient convection-diffusion equations. Applied mathematics 
letters, 18(7):833-841, 2005. 

[46] Alexandre Ern, Annette F. Stephansen, and Martin Vohralik. Guaranteed 
and robust discontinuous Galerkin a posteriori error estimates for 
convection-diffusion-reaction problems. Journal of computational and 
applied mathematics, 234(1):114-130, 2010. 

[47] Don Estep, Michael Pernice, Simon Tavener, and Haiying Wang. A 
posteriori error analysis for a cut cell finite volume method. Computer 
Methods in Applied Mechanics and Engineering, 200(37):2768-2781, 2011. 

[48] L.G. Evans. Partial Differential Eguations. Graduate studies in 
mathematics. American Mathematical Society, 1998. 

[49] Raul Ferreira, Pablo Groisman, and Julio D. Rossi. Numerical blow-up for 
a nonlinear problem with a nonlinear boundary condition. Mathematical 
models and methods in applied sciences, 12(04):461-483, 2002. 

[50] Gadi Fibich and Boaz flan. Discretization effects in the nonlinear 
Schrodinger equation. Applied numerical mathematics, 44(l):63-75, 2003. 

[51] M.H. Friedman. Principles and Models of Biological Transport. Springer, 
2008. 

[52] Victor A. Galaktionov and Juan L. Vazquez. The problem of blow-up in 
nonlinear parabolic equations. Discrete and continuous dynamical systems, 
8(2):399-434, 2002. 

[53] E.H. Georgoulis, E. Hall, and G. Makridakis. An a posteriori error bound 
for discontinuous Galerkin approximations of convection-diffusion problems. 
Submitted for publication, 2014. 

[54] Emmanuil H. Georgoulis, Edward Hall, and Paul Houston. Discontinuous 
Galerkin methods for advection-diffusion-reaction problems on 
anisotropically rehned meshes. SIAM Journal on Scientific Computing, 
30(1):246-271, 2007. 


105 



[55] Emmanuil H. Georgoulis and Omar Lakkis. A posteriori error bounds 
for discontinuous Galerkin methods for quasilinear parabolic problems. In 
Numerical Mathematics and Advanced Applications 2009, pages 351-358. 
Springer, 2010. 

[56] Emmanuil H. Georgoulis, Omar Lakkis, and Juha M. Virtanen. A posteriori 
error control for discontinuous Galerkin methods for parabolic problems. 
SIAM J. Numer. Anal, 49(2):427-458, 2011. 

[57] Emmanuil H. Georgoulis and Gharalambos Makridakis. On a posteriori 
error control for the Allen-Gahn problem. Mathematical Methods in the 
Applied Sciences, 37(2):173-179, 2014. 

[58] Pierre Grisvard. Elliptic problems in nonsmooth domains, volume 69. SIAM, 
2011 . 

[59] J.S. Hesthaven and T. Warburton. On the constants in hp-£nite element 
trace inverse inequalities. Computer methods in applied mechanics and 
engineering, 192(25):2765-2773, 2003. 

[60] Ghiaki Hirota and Kazufumi Ozawa. Numerical method of estimating the 
blow-up time and rate of the solution of ordinary differential equations — an 
application to the blow-up problems of partial differential equations. Journal 
of Computational and Applied Mathematics, 193(2):614 - 637, 2006. 

[61] R. Hoppe, G. Kanschat, and T. Warburton. Gonvergence analysis of an 
adaptive interior penalty discontinuous Galerkin method. SIAM Journal on 
Numerical Analysis, 47(l):534-550, 2009. 

[62] Paul Houston, Ghristoph Schwab, and Endre Siili. Discontinuous hp-hnite 
element methods for advection-diffusion-reaction problems. SIAM Journal 
on Numerical Analysis, 39(6):2133-2163, 2002. 

[63] Paul Houston and Endre Siili. Adaptive Lagrange-Galerkin methods 
for unsteady convection-diffusion problems. Mathematics of computation, 
70(233):77-106, 2001. 


106 



[64] B. Hu. Blow-up Theories for Semilinear Paraholie Equations. Number 
no. 2018 in Blow-up Theories for Semilinear Parabolic Equations. Springer, 
2011 . 

[65] Weizhang Huang, Jingtang Ma, and Robert D. Russell. A study of moving 
mesh PDE methods for numerical simulation of blowup in reaction diffusion 
equations. Journal of Computational Physies, 227(13):6532-6552, 2008. 

[66] B. Janssen and T. P. Wihler. Existence Results for the Continuous and 
Discontinuous Galerkin Time Stepping Methods for Nonlinear Initial Value 
Problems. ArXiv e-prints, July 2014. 

[67] Volker John and Julia Novo. A robust SUPG norm a posteriori error 
estimator for stationary convection-diffusion equations. Computer Methods 
in Applied Meehanies and Engineering, 255:289-305, 2013. 

[68] Claes Johnson, Yi-Yong Nie, and Vidar Thomee. An a posteriori error 
estimate and adaptive timestep control for a backward Euler discretization 
of a parabolic problem. SIAM journal on numerieal analysis, 27(2):277-291, 
1990. 

[69] Chang-Yeol Jung and Roger Temam. Numerical approximation of 
two-dimensional convection-diffusion equations with multiple boundary 
layers. Int. J. Numer. Anal. Model, 2(4):367-408, 2005. 

[70] Ohannes A. Karakashian and Frederic Pascal. A posteriori error 
estimates for a discontinuous Galerkin approximation of second-order 
elliptic problems. SIAM J. Numer. Anal, 41(6):2374-2399 (electronic), 
2003. 

[71] Ohannes A. Karakashian and Frederic Pascal. Convergence of adaptive 
discontinuous Galerkin approximations of second-order elliptic problems. 
SIAM Journal on Numerieal Analysis, 45(2):641-665, 2007. 

[72] Blent Karasozen and Murat Uzunca. Time-space adaptive discontinuous 
Galerkin method for advection-diffusion equations with non-linear reaction 


107 



mechanism. GEM - International Journal on Geomathematics, pages 1-34, 
2014. 


[73] O.T. Kedem and A. Katchalsky. Thermodynamic analysis of the 
permeability of biological membranes to non-electrolytes. Biochimica et 
biophysica Acta, 27:229-246, 1958. 

[74] Daniel Kessler, Ricardo H. Nochetto, and Alfred Schmidt. A posteriori error 
control for the Allen-Cahn problem; circumventing Gronwalls inequality. 
ESAIM: Mathematical Modelling and Numerical Analysis - Modlisation 
Mathmatique et Analyse Numrique, 38(1); 129-142, 2004. 

[75] Christian Klein, Benson Muite, and Kristelle Roidot. Numerical study of 
blowup in the Davey-Stewartson system. arXiv preprint arXiv:l 112.4043, 
2011 . 

[76] Natalia Kopteva and Eugene O’Riordan. Shishkin meshes in the numerical 
solution of singularly perturbed differential equations. Int. J. Numer. Anal. 
Model, 7(3):393-415, 2010. 

[77] Gerd Kunert. A posteriori error estimation for convection dominated 
problems on anisotropic meshes. Mathematical methods in the applied 
sciences, 26(7):589-617, 2003. 

[78] I. Kyza. A posteriori error estimates for approximations of semilinear 
parabolic and Schrodinqer-type equations. PhD thesis, PhD Thesis, 
University of Crete, 2009. 

[79] Irene Kyza and Charalambos Makridakis. Analysis for time discrete 
approximations of blow-up solutions of semilinear parabolic equations. 
SIAM Journal on Numerical Analysis, 49(l):405-426, 2011. 

[80] P. Lesaint and P.A. Raviart. On a Einite Element Method for Solving the 
Neutron Transport Equation. Univ. Paris VI, Labo. Analyse Numerique, 
1974. 


108 



[81] Charalambos Makridakis and Ricardo H. Nochetto. Elliptic reconstruction 
and a posteriori error estimates for parabolic problems. SIAM J. Numer. 
Anal, 41 (4): 1585-1594, 2003. 

[82] P. Morin, R. Nochetto, and K. Siebert. Data oscillation and convergence of 
adaptive FEM. SIAM Journal on Numerical Analysis, 38(2):466-488, 2000. 

[83] Lin Mu and Rabeea Jari. A posteriori error analysis for discontinuous finite 
volume methods of elliptic interface problems. Journal of Computational 
and Applied Mathematics, 255:529-543, 2014. 

[84] Gohisse N., Firmin K., and Theodore K. Boni. Numerical blow-up for 
a nonlinear heat equation. Acta Mathematica Sinica, English Series, 
27(5):845-862, 2011. 

[85] V.T. Nguyen and H. Zaag. Blow-up results for a strongly perturbed 
semilinear heat equation: Theoretical analysis and numerical method. 2014. 

[86] Louis Nirenberg. On elliptic partial differential equations. Springer, 2011. 

[87] J. Nitsche. Uber ein variationsprinzip zur losung von dirichlet-problemen 
bei verwendung von teilraumen, die keinen randbedingungen unterworfen 
sind. In Ahhandlungen aus dem mathematischen Seminar der Universitdt 
Hamburg, volume 36, pages 9-15. Springer, 1971. 

[88] M. Picasso. Adaptive hnite elements for a linear parabolic 
problem. Computer Methods in Applied Mechanics and Engineering, 
167(3-4) :223-237, 1998. 

[89] Marco Picasso and Virabouth Prachittham. An adaptive algorithm for the 

Crank-Nicolson scheme applied to a time-dependent convection-diffusion 
problem. Journal of computational and applied mathematics, 

233(4):1139-1154, 2009. 

[90] Dirk Praetorius, Ewa Weinmiiller, and Philipp Wissgott. A space-time 
adaptive algorithm for linear parabolic problems. 


109 



[91] Martin Prosi, Paolo Zunino, K. Perktold, and Alfio Quarteroni. 
Mathematical and numerical models for transfer of low-density lipoproteins 
through the arterial walls: a new methodology for the model set up 
with applications to the study of disturbed lumenal flow. Journal of 
biomechanics, 38(4):903-917, 2005. 

[92] W.H. Reed and T.R. Hill. Triangular mesh methods for the neutron 
transport equation. Los Alamos Report LA-UR-73-479, 1973. 

[93] H.G. Roos, M. Stynes, and L. Tobiska. Robust Numerical Methods for 
Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction 
and Flow Problems. Springer Series in Computational Mathematics. 
Springer, 2008. 

[94] Giancarlo Sangalli. A uniform analysis of nonsymmetric and coercive linear 
operators. SIAM journal on mathematical analysis, 36(6):2033-2048, 2005. 

[95] Giancarlo Sangalli. Robust a-posteriori estimator for 

advection-diffusion-reaction problems. Mathematics of Computation, 
77(261):41-70, 2008. 

[96] Alfred Schmidt and Kunibert G. Siebert. Design of adaptive finite 
element software, volume 42 of Lecture Notes in Computational Science 
and Engineering. Springer-Verlag, Berlin, 2005. The hnite element toolbox 
ALBERTA, With 1 CD-ROM (Unix/Linux). 

[97] Dominik Schotzau and Liang Zhu. A robust a-posteriori error estimator for 
discontinuous Galerkin methods for convection-diffusion equations. Appl. 
Numer. Math., 59(9):2236-2255, 2009. 

[98] A.M. Stuart and M.S. Floater. On the computation of blow-up. Euro. J. 
Appl. Math., 1:47-71, 1990. 

[99] Shuyu Sun and Mary F. Wheeler. A posteriori error estimation and 
dynamic adaptivity for symmetric discontinuous Galerkin approximations 
of reactive transport problems. Computer methods in applied mechanics and 
engineering, 195(7):632-652, 2006. 


110 



[100] T. Teorell. Transport processes and electrical phenomena in ionic 
membranes. Prog. Biophys. Biophys. Chem., 3:305-369, 1953. 

[101] Y. Tourigny and J.M. Sanz-Serna. The numerical study of blowup with 
application to a nonlinear Schrodinger equation. Journal of Computational 
Physics, 102(2) :407-416, 1992. 

[102] Takeo K. Ushijima. On the approximation of blow-up time for solutions 
of nonlinear parabolic equations. Publications of the Research Institute for 
Mathematical Sciences, 36(5):613-640, 2000. 

[103] R. Verfiirth. A posteriori error estimates for non-linear parabolic equations. 
Preprint, Ruhr-Universitdt Bochum, Fakultdt fiir Mathematik, Bochum, 
Germany, 2004. 

[104] R. Verfiirth. Robust a posteriori error estimates for nonstationary 
convection-diffusion equations. SIAM J. Numer. Anal, 43(4): 1783-1802 
(electronic), 2005. 

[105] Riidiger Verfiirth. A Review of A Posteriori Error Estimation and Adaptive 
Mesh-Refinement Techniques. Wiley-Teubner, Chichester-Stuttgart, 1996. 

[106] Riidiger Verfiirth. A posteriori error estimates for nonlinear problems: 
L''(0, f; L^(r2))-error estimates for hnite element discretizations of parabolic 
equations. Mathematics of Computation of the American Mathematical 
Society, 67(224):1335-1360, 1998. 

[107] Riidiger Verfiirth. A posteriori error estimates for nonlinear problems: 
L^(0, t; lV^’^(r2))-error estimates for finite element discretizations of 
parabolic equations. Numerical Methods for Partial Differential Equations, 
14(4):487-518, 1998. 

[108] Riidiger Verfiirth. A posteriori error estimators for convection-diffusion 
equations. Numerische Mathematik, 80(4):641-663, 1998. 

[109] Riidiger Verfiirth. Error estimates for some quasi-interpolation operators. 
ESAIM: Mathematical Modelling and Numerical Analysis, 33(04):695-713, 
1999. 


Ill 



[110] Riidiger Verfiirth. Robust a posteriori error estimates for stationary 
convection-diffusion equations. SIAM journal on numerical analysis, 
43(4): 1766-1782, 2005. 

[111] Mary Fanett Wheeler. A priori error estimates for Galerkin 
approximations to parabolic partial differential equations. SIAM Journal 
on Numerical Analysis, 10(4):723-759, 1973. 

[112] Weiying Zheng and He Qi. On Friedrichs-Poincare-type inequalities. 
Journal of mathematical analysis and applications, 304(2):542-551, 2005. 

[113] Liang Zhu and Dominik Schotzau. A robust a posteriori error estimate for 
hp-adaptive dG methods for convection-diffusion equations. IMA journal of 
numerical analysis, 31(3):971-1005, 2011. 

[114] Paolo Zunino. Mathematical and numerical modeling of mass transfer in 
the vascular system. PhD thesis, Politecnico di Milano, 2002. 


112 



